Atmospheric correction algorithm for short-wave channels of the MSU-MR device of the Meteor-M no. 2 satellite

Full Text

ВВЕДЕНИЕ

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

АК спутниковых измерений является необходимым условием успешного решения широкого спектра гидрологических и океанографических задач, например, расчета вегетационных индексов (Додышева, 2010; Зубкова и др., 2015), поскольку даже незамутненная атмосфера оказывает значительное влияние на их значение. В Дальневосточном центре НИЦ «Планета» при получении характеристик подстилающей поверхности длительное время использовались данные зарубежных спутниковых систем, которые предоставляются ведущими мировыми космическими агентствами уже скорректированными. В настоящее время для получения продуктов с российских полярно-орбитальных космических аппаратов (КА) используются данные с двух спутников («Метеор-М» № 2 и «Метеор-М» № 2–2), а до 2024 года предполагается наличие на орбите уже шести КА этой серии. Основная проблема в настоящее время заключается в отсутствии системного подхода к проведению АК для данных многозонального сканирующего устройства малого разрешения (МСУ-МР), установленного на борту спутников этой серии. Такая ситуация не позволяет полноценно использовать эти данные в задачах гидрометеорологии в виду несоответствия получаемых параметров требованиям Всемирной метеорологической организации (ВМО), предъявляемым к точности спутниковых продуктов.

Целью настоящей работы является реализация программных средств АК спутниковых измерений в коротковолновой части спектра, которые позволят учитывать максимальное количество искажающих факторов. Большинство существующих программных средств АК основываются на справочных таблицах (Look-Up Table, LUT), содержащих смоделированные значения интенсивности излучения в каналах спутникового прибора для различной геометрии наблюдения и состояния атмосферы. Данные таблицы строятся с помощью моделей переноса излучения (например, MODTRAN или 6S). По результатам сравнительного анализа для решения задачи АК измерений коротковолновых каналов МСУ-МР выбрана модель переноса излучения Second Simulation of a Satellite Signal in the Solar Spectrum (6S) (Kotchenova, Vermote, 2007). Разработана функция для генерирования LUT с использованием выбранной модели. Реализован алгоритм восстановления скорректированных значений КСЯ с использованием LUT. Проведена верификация первого канала МСУ-МР, показывающая высокую корреляцию с эталонными значениями, в качестве которых выступают КСЯ подстилающей поверхности тестовых полигонов по данным портала EUMETSAT «Surface Albedo Validation Sites» (EUMETSAT, 2015). Дополнительная проверка реализованного алгоритма применительно к данным первого канала прибора Advanced Very High Resolution Radiometer (AVHRR) КА «MetOp-A» показала, что корреляция КСЯ с эталонным значением для обоих приборов сопоставима.

ОБЗОР ИССЛЕДОВАНИЙ И РАБОТ

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

Физические алгоритмы АК в большинстве своем основаны на предварительно рассчитанных LUT, которые строятся с использованием модели переноса излучения. LUT представляют собой таблицы с заранее смоделированными выходными значениями и параметрами, которые закладывались при моделировании. Такой подход значительно сокращает время работы алгоритмов, поскольку поиск ближайших значений в таблице происходит намного быстрее самого процесса моделирования. В настоящее время LUT таблицы для различных алгоритмов моделируются, в основном, либо с использованием модели переноса излучения Moderate resolution atmospheric transmission (MODTRAN) (Berk et al., 2014), либо с использованием модели 6S (Kotchenova, Vermote, 2007). Обе модели схожи по входным параметрам, отличие заключается в самом алгоритме моделирования и режимах работы. 6S позволяет работать в двух режимах: имитация, когда модель генерирует атмосферный фон и значение КСЯ ВГА со всеми остальными параметрами для решения задачи АК, и корректирование, при котором модель, учитывая КСЯ ВГА и параметры атмосферы, вычисляет КСЯ подстилающей поверхности. MODTRAN работает только в режиме корректирования.

В свою очередь, математические алгоритмы АК представлены как простыми реализациями, например, Dark-Object Subtraction (DOS) (Liang et al., 2001), так и усложненными, к примеру, Quick Atmospheric Correction (QUAC) (Bernstein et al., 2014). С использованием математических алгоритмов при проведении АК нет необходимости задавать оптические параметры аэрозоля, общее содержание газов и учитывать подстилающую поверхность, не требуется наличие предварительно рассчитанных LUT. Однако применение этих алгоритмов имеет ряд ограничений. К примеру, DOS, основанный на поиске минимального значения в области интереса с последующим его вычитанием из значения каждого пикселя, требует присутствия в сцене темного объекта (обычно это вода), а также применяется только на малых участках. QUAC основан на эмпирическом выводе о том, что средняя отражательная способность спектров различных поверхностей не зависит от каждой сцены, однако это не всегда соответствует действительности и может вносить погрешность в АК. Тем не менее, в настоящее время такие алгоритмы активно развиваются и их точность может быть соизмерима с физическими. В работе (Беляев и др., 2016) разработан алгоритм АК применительно к данным ручных видеоспектральных приборов, установленных на борту российского сегмента Международной космической станции (МКС). Так же предложен алгоритм коррекции влияния газового поглощения на КСЯ, который не требует решения сложной оптимизационной задачи (Николаева, 2016). Особенностями обоих алгоритмов является использование значений в узких спектральных интервалах (от 5 до 22 нм), что с одной стороны позволяет добиться высокой точности, а с другой накладывает ограничения на их применимость.

В работе (Nguyen et al., 2015) определялся оптимальный метод АК для оценки лесной биомассы по данным прибора Enhanced Thematic Mapper Plus (ETM+) КА «Landsat». Для этого тестировались три наиболее популярных алгоритма: Fast Line-of-sight Atmospheric Analysis of Hypercubes (FLAASH) (Adler-Golden et al., 1998), основанный на LUT по данным модели MODTRAN и требующий ручного ввода характеристик изображения и знания параметров атмосферы; DOS; модель 6S. Была проведена серия экспериментов для оценок лесной биомассы для территории Южной Кореи, в результате которых модель 6S показала наименьшие значения среднеквадратической ошибки (RMSE). В дополнение, по итогам проверки чувствительности методов АК к каждой из полос ETM+, 6S показала лучшие результаты.

Целью проекта (Smith, 2015) являлось сравнение трех алгоритмов АК применительно к спутниковым данным КА «WorldView‑2»: DigitalGlobe Atmospheric Compensation (DG AComp) (Pacifici et al., 2014), основанный на данных модели MODTRAN и автоматически определяющий атмосферные условия; QUAC, который использует математический подход; FLAASH. Результаты оценки точности в этом исследовании свидетельствуют о том, что DG AComp является наиболее корректным алгоритмом АК, он почти в два раза точнее, чем FLAASH. Кроме того, DG AComp обладает наилучшей точностью для каждого типа поверхности и спектрального диапазона. QUAC, в свою очередь, имеет наихудшую точность из всех методов.

Результаты (Kotchenova et al., 2008) представлены в научном проекте, посвященном сравнению четырех моделей переноса излучения в атмосфере: 6S, MODTRAN, Spherical Harmonics (SHARM) (Muldashev et al., 1999) и Radiative Transfer (RT3) (Evans, Stephens, 1991). Надежность моделей проверялась с помощью известных эталонных тестов, таких как код Монте Карло (Monte Carlo) (Breon, 2014) и табличные значения Коулсона (Coulson) (Coulson et al., 1961). Сравнение проводилось отдельно для молекулярной, аэрозольной и смешанной атмосферы. 6S продемонстрировала лучшее согласование с Coulson и Monte Carlo. Более того, при дальнейших сравнениях 6S использовалась в качестве эталона вместо Monte Carlo, поскольку она способна обеспечивать моделирование, которое является столь же точным.

Ввиду того, что точность математических алгоритмов в рассмотренных исследованиях (Nguyen et al., 2015; Smith, 2015; Kotchenova et al., 2008) проигрывает физическим, последние было решено использовать для задачи АК применительно к данным МСУ-МР. Поскольку модель 6S в (Nguyen et al., 2015; Kotchenova et al., 2008) имеет наибольшую точность, а в (Kotchenova et al., 2008) использовалась в качестве эталона, она была выбрана в настоящей работе для моделирования переноса излучения.

ГЕНЕРИРОВАНИЕ И ИСПОЛЬЗОВАНИЕ LUT

Моделирование с использованием 6S ограничено длинами волн от 0.4 до 4.0 мкм, поэтому в работе рассматривались функции спектральной чувствительности только коротковолновых каналов прибора МСУ-МР, показанных на рис. 1.

 

Рис. 1. Функции спектральной чувствительности коротковолновых каналов МСУ-МР.

 

Для решения задачи АК коротковолновых каналов прибора МСУ-МР с использованием модели 6S необходимо, в первую очередь, сгенерировать LUT. Для успешного моделирования необходимо задать ряд параметров, которые заранее определены в самой модели. Для этих параметров в определенном диапазоне составляются дискретные наборы значений, запускается моделирование с использованием 6S для каждой возможной комбинации, после чего результат и соответствующие входные параметры сохраняются в таблице.

Функции спектральной чувствительности каналов прибора МСУ-МР задаются с шагом 0.025 мкм и, также как стартовая и конечная длины волн каналов, являются неизменными в процессе моделирования. Таблицы для разных каналов моделируются отдельно друг от друга.

Для генерирования LUT была написана функция, с использованием библиотеки Py6S (Wilson, 2013), в которой задаются параметры, определенные в 6S:

— КСЯ подстилающей поверхности (от 0 до 1 с шагом 0.0125);

— солнечный зенитный угол (от 0° до 80° с шагом 10°);

— зенитный угол наблюдения (от 0° до 60° с шагом 10°);

— относительный азимутальный угол (от 0° до 180° с шагом 60°);

— высота поверхности относительно уровня моря (от 0 км до 9 км с шагом 3 км);

— оптическая толщина аэрозоля (AOT) на длине волны 550 нм (от 0 до 0.3 с шагом 0.05, от 0.3 до 5 с шагом 0.2);

— содержание водяного пара в вертикальном столбе атмосферы (от 0 г/см² до 10 г/см² с шагом 1 г/см²);

— содержание озона в вертикальном столбе атмосферы (от 100 DU (Dobson unit) до 500 DU с шагом 100 DU);

— аэрозольный профиль, который предустановлен в модели 6S в значениях «континентальный» или «морской».

Стоит отметить, что выбранный шаг моделирования параметров обусловлен временными затратами на генерирование таблицы. Так, при генерировании LUT для одного канала МСУ-МР с использованием описанной функции применялось 4 сервера: два с процессорами Intel Xeon E5–2687W v4 и два — с Intel Xeon Gold 6130, при этом время на генерирование в этом случае составило 4 месяца.

Результатом выполнения функции является получение значения эффективной СПЭЯ ВГА для каждого канала прибора МСУ-МР в отдельности. После проведения моделирования набор входных и выходных параметров заносится в базу данных для быстрого поиска записей.

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

При АК для работы с LUT применяются СПЭЯ, рассчитываемые по формуле (1) с использованием значений КСЯ, прошедших интеркалибровку по данным прибора AVHRR (Филей и др., 2016).

$СПЭ{Я}_i=\frac{\left(КС{Я}_i·F_i·\cos \left(\Theta \right)\right)}{\left(\pi ·d^2·100\right)}$ , (1)

где СПЭЯ_i — СПЭЯ i-го канала, КСЯ_i — КСЯ (в процентах) i-го канала, F_i — интегральный поток солнечного излучения, свернутый с относительной спектральной функцией прибора, i — порядковый номер канала прибора МСУ-МР,  — зенитный угол Солнца, d — среднее расстояние от Земли до Солнца в астрономических единицах.

Данные об углах Солнца и наблюдения рассчитываются по имеющейся информации о положении КА «Метеор-М» № 2 в момент съемки. Данные о содержании водяного пара и озона берутся из прогностической модели GFS (Global Forecast System). Для совмещения времени съемки и данных модели применяется интерполяция по времени с использованием оптического потока (Симоненко и др., 2017). В качестве аэрозольного профиля для суши используется континентальный профиль, для воды — морской. Высота поверхности над уровнем моря берется из ресурса NASA Shuttle Radar Topographic Mission (SRTM) (Jarvis et al., 2008). АК проводится только для безоблачных пикселей, определенных по маске облачности с применением сверточной нейронной сети (Андреев и др., 2019). Поскольку в распоряжении авторов отсутствует технология восстановления АОТ по данным прибора МСУ-МР, то при проведении АК используются данные прибора MODIS КА «TERRA» и «AQUA», продукт MOD04 (Aerosol Product) (Kaufman et al., 1997). Вследствие того, что используемые данные имеют разное пространственное разрешение, они приводятся к разрешению каналов прибора МСУ-МР, которое составляет 1 км.

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

Алгоритм для группы пикселей находит ближайшие значения всех параметров из LUT и считывает значения смоделированных СПЭЯ в паре с КСЯ. После этого для каждой группы производится сначала интерполяция СПЭЯ по реальным данным с использованием кубических сплайнов (Самарский, Гулин, 1989), а затем процесс повторяется итерационно для каждого параметра с использованием линейной интерполяции (2):

$Alb=Al{b}_1+\left(Al{b}_2-Al{b}_1\right)·\left(\frac{Va{l}_{curr}-Va{l}_1}{Va{l}_2-Va{l}_1}\right)$ (2)

где Alb — результирующее значение КСЯ, Val₁ и Val₂ — граничные значения параметра, Alb₁ и Alb₂ — значения КСЯ для Val₁ и Val₂, Val_curr — текущее значение параметра в пикселе.

ВЕРИФИКАЦИЯ РЕЗУЛЬТАТОВ АК

Верификация результатов АК производилась с использованием базы данных портала EUMETSAT «Surface Albedo Validation Sites» (EUMETSAT, 2015). На портале данные о КСЯ поверхности в видимом спектре предоставлены для конкретной области (полигона) по информации с одной из следующих сетей: FLUXNET, ALBEDOVAL, BELMANIP, AERONET и др. В работе КСЯ поверхности таких областей использовалось как эталонное. Из базы данных было отобрано 30 однородных участков земной поверхности, КСЯ которых не подвергалось годовой и сезонной изменчивости за весь период наблюдения. На рис. 2 показаны некоторые из полигонов по данным сервиса GoogleMaps. Ввиду того, что значения КСЯ для полигонов предоставлены только в видимом спектре, верификация проводилась только для первого канала прибора МСУ-МР.

Для дополнительной проверки посредством разработанного алгоритма атмосферная коррекция была проведена для данных первого канала прибора AVHRR КА «MetOp-A», имеющего схожие спектральные характеристики с первым каналом прибора МСУ-МР (рис. 3).

Для верификации были подобраны пролеты КА «Метеор-М» № 2 и «MetOp-A» в те дни, когда над тестовыми полигонами отсутствовала облачность. Для первого канала прибор AVHRR, с использованием его функции спектральной чувствительности (рис. 3), начальной и конечной длин волн, была сгенерирована отдельная LUT через Py6S. При генерации данной LUT диапазон каждого параметра был ограничен его максимальным и минимальным значением на тестовых полигонах, что позволило уменьшить временные затраты при моделировании.

 

Рис. 2. Тестовые полигоны.

  

Рис. 3. Функции спектральной чувствительности первых каналов приборов AVHRR и МСУ-МР.

 

 

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

$КС{Я}_{ср}=\left(\frac{1}{n}\right)·\sum _{i=1}^nКС{Я}_i$ , (3)

где КСЯ_ср — результирующее среднее значение КСЯ, КСЯ_i — значение КСЯ в i-ом пикселе тестового полигона, n — количество пикселей тестового полигона.

Для значений КСЯ_ср, которые выражены в процентах, в табл. 1 приведены абсолютные ошибки, а на рис. 4 они соотнесены с эталонным КСЯ полигонов (КСЯ^ref). Средняя ошибка (E_mean) и среднеквадратичная ошибка (E_RMS) в табл. 1 рассчитываются по формулам (4) и (5) соответственно. Коэффициент детерминации (R) рассчитывается согласно (Бахрушин, 2011).

$E_{mean}=\left(\frac{1}{n}\right)·\sum _{i=1}^n\left(КС{Я}_i^{ref}-КС{Я}_i\right)$ , (4)

где E_mean — средняя ошибка, – эталонное значение КСЯ i-го тестового полигона, КСЯ_i — значение КСЯ_ср i-го тестового полигона, n — количество тестовых полигонов

$E_{RMS}=\sqrt{\left(\frac{1}{n}\right)·\sum _{i=1}^n{\left(КС{Я}_i^{ref}-КС{Я}_i\right)}^2}$ , (5)

где E_RMS — среднеквадратичная ошибка, $КС{Я}_i^{ref}$ – эталонное значение КСЯ i-го тестового полигона, КСЯ_i — значение КСЯ_ср i-го тестового полигона, n — количество тестовых полигонов.

Анализируя E_mean и E_RMS из табл. 1 и значения R из рис. 4 можно сделать вывод, что для обоих спутниковых приборов после проведения АК значения КСЯ_ср стали в меньшей степени отклоняться от КСЯ^ref. Исходя из вышесказанного, можно сделать вывод о том, что АК сопоставимо увеличивает корреляцию с эталонными значениями КСЯ для обоих спутниковых приборов.

 

Таблица 1. Абсолютные ошибки для КСЯ, выраженных в процентах

Абсолютные ошибки	МСУ-МР		AVHRR
	до коррекции	после коррекции	по коррекции	после коррекции
Emean	1.594	0.468	2.587	1.074
ERMS	2.024	1.098	1.967	1.556
R	0.88	0.93	0.819	0.923

 

Рис. 4. Соотношение КСЯ до АК и после нее с эталонным КСЯ (1 — для первого канала прибора МСУ-МР, 2 — для первого канала прибора AVHRR).

ЗАКЛЮЧЕНИЕ

В результате работы был реализован алгоритм АК с использованием модели переноса излучения 6S применительно к данным прибора МСУ-МР. По результатам верификации первого канала прибора, при проведении которой были рассчитаны средние и среднеквадратичные ошибки, а также коэффициент детерминации, можно сделать вывод, что АК увеличивает корреляцию значений КСЯ по данным МСУ-МР с эталонными, в качестве которых использовались КСЯ подстилающей поверхности тестовых полигонов по данным портала EUMETSAT «Surface Albedo Validation Sites». При сравнении результатов работы алгоритма для данных прибора AVHRR и МСУ-МР корреляция в обоих случаях сопоставима.

При использовании алгоритма в оперативной работе есть ограничение, связанное с использованием данных АОТ прибора MODIS КА «TERRA» и «AQUA». Для проведения АК в этом случае необходимо либо искать пересекающие пролеты спутников, либо использовать глобальные суточные карты АОТ. Для решения этой проблемы в настоящий момент в Дальневосточном центре НИЦ «Планета» ведутся работы по восстановлению параметра АОТ по данным прибора МСУ-МР.

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

В дальнейших планах стоит задача реализации верификации настоящего алгоритма АК для второго и третьего каналов прибора МСУ-МР, а также реализация алгоритма проведения АК с использованием нейросетевого подхода.

БЛАГОДАРНОСТИ

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

About the authors

M. O. Kuchma

Far-Eastern Center of State Research Center for Space Hydrometeorology «Planeta»

Author for correspondence.
Email: m.kuchma@dvrcpod.ru
Russian Federation, Khabarovsk

V. D. Bloshchinskiy

Far-Eastern Center of State Research Center for Space Hydrometeorology «Planeta»

Email: m.kuchma@dvrcpod.ru
Russian Federation, Khabarovsk

References

Supplementary files