Развитие методов предварительной обработки данных MODIS, VIIRS
Данные дистанционного зондирования Земли позволяют осуществлять непрерывный глобальный мониторинг состояния растительного покрова. Для анализа вегетационных характеристик распространено использование наблюдений, полученных в красном (RED), ближнем инфракрасном (NIR), коротковолновом инфракрасном (SWIR) спектральных диапазонах, и построенных на их основе индексах. На информативность большинства данных, отражающих состояние растительного покрова влияют такие факторы, как облачный и снежный покров, тени от облаков, атмосферная дымка и т.д.
В ИКИ РАН разработана и используется технология предварительной обработки данных спутниковых наблюдений, полученных приборами Terra/Aqua-MODIS и NPP-VIIRS. Предварительная обработка данных включает в себя формирование маскирующих изображений облачного и снежного покрова, построение теней от облаков, фильтрацию наблюдений, полученных при некорректных углах съёмки и создание разновременных композитных изображений. В результате фильтрации измерений, искажённых под влиянием мешающих факторов, в изображениях остаются пропуски, остаточные шумы и выбросные значения. Разработанный в ИКИ РАН интерполяционный алгоритм восстановления длинных временных рядов данных спутникового наблюдения позволяет восстанавливать пропуски, корректировать остаточные искажения и сглаживать временные ряды наблюдений.
На первом этапе предварительной обработки происходит формирование маскирующих изображений облачного и снежного покрова, включающих в себя все неудовлетворительные измерения сеансов наблюдений. Для создания масок облачного и снежного покрова используются значения коэффициентов спектральной яркости (КСЯ) в каналах b01, b02, b03 и b06 для MODIS и I-1, I-2, M-3, I-3 для VIIRS.
Таблица 1. Спектральные каналы VIIRS и MODIS для построения масок облачного и снежного покрова
|
|
MODIS
|
VIIRS
|
|
band
|
λ, мкм
|
band
|
λ, мкм
|
|
RED
|
b01
|
0,620–0,670
|
I-1
|
0,597–0,679
|
|
NIR
|
b02
|
0,841–0,876
|
I-2
|
0,842–0,881
|
|
BLUE
|
b03
|
0,459–0,479
|
M-3
|
0,477–0,496
|
|
SWIR
|
b06
|
1,628–1,652
|
I-3
|
1,570–1,629
|
При детектировании снежного и облачного покрова используется разностный нормализованный индекс снега NDSI. Для расчета NDSI могут браться разные спектральные каналы. Как правило индекс отражает разность КСЯ в одном из видимых спектральных каналов и КСЯ коротковолнового инфракрасного диапазона. При создании маски облачного покрова применяются индексы NDSI, построенные с использованием красного по формуле (1) и голубого (2) каналов:
(1)
(2)
Значения индексов NDSI(R) и NDSI(B), а также КСЯ в каналах, перечисленных в таблице 1, используются для предварительного порогового маскирования облачного и снежного покрова.
На рисунке 2 приведён пример построения предварительной маски облачного и снежного покрова. Пример содержит множество мелких облаков, размер которых сопоставим с пространственным разрешением данных. Такие облака сложно детектировать из-за спектрального смешения с подстилающей поверхностью. Для фильтрации смешанных пикселей на предварительном этапе используется оконтуривание матрицей 3×3 пикселя для облаков типа CLD и SLD (рисунок 3).
Таблица 2. Набор пороговых значений для предварительного маскирования облаков
|
Маркер
|
Характеристика
|
Условия маркирования
|
|
BAD
|
Темные объекты: водная поверхность, тени, ошибки
|
B01 < 0, B02 < 0
B03 < 0, B06 < 0
|
|
SNOW
|
Снег, сильная облачность
|
B01 > 700, NDSI(R) > 0,1
B03 > 700, NDSI(B) > 0,1
|
|
CLD
|
Высокая облачность
|
B01 > 700, NDSI(R) > –0,2
B03 > 700, NDSI(B) > –0,2
Оконтуривание: 3×3 пикселя
|
|
SLD
|
Средняя облачность
|
B01 > 700, NDSI(R) > –0,35
B03 > 700, NDSI(B) > –0,35
Оконтуривание: 3×3 пикселя
|
|
MLD
|
Дымка, смешанные пиксели
|
B03 > 700, NDSI(B) > –0,45
|
Рисунок 1 — Исходные данные. RGB синтез каналов: R — b02, G — b06, B 3 — b01
Рисунок 2 — Предварительное маскирование облачного покрова. Тёмно-синий цвет — облака типа CLD, синий — облака типа SLD, голубой — облака типа MLD
Рисунок 3 — Оконтуривание облаков типа CLD и SLD. Голубой цвет — все маркеры предварительной маски облачного покрова, синий — оконтуривание облаков
Рисунок 4 — Предварительное маскирование облачного покрова
Недостатком порогового метода является высокий процент ошибок детектирования: ложное детектирование и пропуск цели. Возможность достоверного разделения чистых наблюдений и облачного покрова варьируется в зависимости от специфики и плотности облаков, времени года, типа и состояния растительного покрова. Для уточнения маскирующих изображений и повышения адаптивности метода в целом используется пространственный гистограммный анализ значений КСЯ в голубом канале и индекса NDSI.
Для каждого класса облаков строится буферная зона (рисунок 5). Гистограммы распределения значений внутри класса «облако» (рисунок 7) и «буфер» (рисунок 8) частично накладываются. На рисунке 9 представлена нормализованная разница гистограмм распределения значений КСЯ в голубом канале для облаков типа SLD и соответствующего буфера.
Рисунок 5 – Построение буферной зоны (оранжевый) вокруг облаков типа CLD(синий)
Рисунок 6 – Коррекция (красный) облачного покрова
Рисунок 7 — Гистограмма распределения значений в голубом канале b03 для облаков типа SLD
Рисунок 8 — Гистограмма распределения значений в голубом канале b03 для буферной зоны
Рисунок 9 — Разностная гистограмма объектов класса «облако» и «буфер»
Точка перемены знака разностной гистограммы позволяет уточнить порог разделения классов для конкретного изображения. Полученный порог рассчитывается дважды: для фильтрации ошибок ложного детектирования внутри класса «облако», и, после перестроения буферной зоны, для детектирования облачности внутри буфера.
На рисунке 6 приведена маска облачного покрова после пространственной гистограммной коррекции. Красным цветом на ней выделены облака, которые удалось добавить к предварительной маске облачного покрова.
Следующий этап предварительной обработки данных спутниковых наблюдений включает в себя фильтрацию наблюдений, полученных при величине зенитного угла Солнца выше 75° и величине зенитного угла наблюдения выше 44° для данных MODIS и выше 67° для данных VIIRS, а также построение теней от облаков. При этом также используется информация о геометрии наблюдения. Направление теней вычисляется исходя из положения Солнца и спутника во время съемки. Длина участка, возможного для присутствия теней, зависит от высоты облачного покрова. На предварительном этапе при построении высота облаков принимается равной 8000 м.
На рисунке 10 приведён пример построения теней с использованием маски облачного покрова и данных об азимутальных и зенитных углах Солнца и сенсора. При этом область потенциальных теней заведомо избыточна и может превышать реальную в несколько раз.
Приведённый выше пространственный гистограммный анализ используется для коррекции класса «тени». Для класса «тени» также строится буферная зона (рисунок 11). Разностные гистограммы распределения значений КСЯ в NIR и SWIR каналах позволяют получить пороги разделения классов «тени» и «подстилающая поверхность» для выбранного изображения. При этом исправляются только ошибки ложного детектирования.
Рисунок 10 — Построение теней (зелёный цвет) от облаков (голубой)
Рисунок 11 — Построение буферной зоны (оранжевый цвет) вокруг теней (зелёный)
Рисунок 12 — Маскирующие изображение облачного и снежного покрова, теней от облаков, шумов и т.д.
Рисунок 13 — Фильтрация искажённых данных (чёрный цвет)
На рисунке 12 приведён окончательный результат построение маскирующего изображения облачного и снежного покрова, теней от облаков, шумов и т.д. Маски используются для фильтрации наблюдений, непригодных для анализа растительного покрова. При этом как видно на рисунке 13, получаются пропуски в данных.
Третьим этапом предварительной обработки является построение композитных изображений. На каждую дату может быть несколько сеансов наблюдений, и создание композитных изображений позволяет частично заполнить «пропуски» в данных. Для различных задач дистанционного зондирования могут использоваться однодневные, четырёхдневные, недельные и месячные композиты.
Для построения композитов используются данные сеансов наблюдения внутри выбранного временного интервала построения композита, а также наблюдения в расширенной временном диапазоне для сбора статистики в каждой точке наблюдения.
На рисунке 14 приведены исходные сеансы наблюдения в RGB синтезе ближнего инфракрасного, коротковолнового и красного каналов и построенные для них маскирующие изображения облачного покрова для ближайших дат к наблюдению на рисунке 13. В двухнедельный промежуток времени один и тот же объект может несколько раз классифицироваться как облако, тень или чистая земная поверхность. Статистический анализ позволяет корректировать классы «облака» и «тени» в масках облачного покрова с учётом значений КСЯ в спектральных каналах.
Для объектов класса «тени» за разные дни вычисляется среднее и стандартное отклонение. Каждое наблюдение, маркированное как «тени», сравнивается со средним значением КСЯ теней в этой точки. При выходе за пределы одного стандартного отклонение наблюдение маркируется как «открытая поверхность». Аналогично проверяются наблюдения открытой поверхности на принадлежность к классу облаков.
На основе чистых наблюдений строится композитное изображение. Если чистых наблюдений больше одного, вычисляется среднее значение и выбирается реальное наблюдение, наиболее близкое к среднему.
Построенные композитные изображения в зависимости от выбранного временного разрешения все ещё содержат пропуски в данных и остаточные искажения. На заключительном этапе предварительной обработки данных происходит анализ временных рядов для фильтрации выбросов, заполнения пропусков и сглаживания данных.
Алгоритм восстановления данных включает в себя построение интерполирующих полиномов второй степени в скользящем окне динамического размера в p значащих точек. Число значащих точек зависит от временного разрешения исходного набора данных, а длина полинома по оси абсцисс зависит от числа пропусков. На рисунке 15 приведён пример в виде набора аппроксимирующих полиномов для временного ряда недельных композитных изображений вегетационного индекса NDVI за один календарный год.
194-a 194-b 201-a 201-b
195-a 195-b 205-a 205-b
199-a 199-b 207-a 207-b
Рисунок 14 — Несколько наблюдений в диапазоне 194–207 день года; a — RGB синтез каналов: R – b02, G — b06, B – b01; b — маски облачного покрова
Рисунок 15 — Интерполяционный алгоритм восстановления длинных временных рядов данных спутникового наблюдения
Результаты интерполяции записываются в интерполяционную матрицу IM. Для каждой точки временного ряда получается набор интерполяционных оценок значения вегетационной характеристики в каждой временной точке: при p = 5 до четырёх оценок для пропуска и до пяти для значащей точки во временном ряду.
(3)
Реальные наблюдения сравниваются с интерполяционными оценками. При отклонении в две дисперсии от среднего наблюдение считается искаженным и фильтруется. После фильтрации выбросов интерполяционная матрица оценок перестраивается, и пропуски заполняются усреднёнными значениями интерполяционных оценок.
В матрице IM элементы aijT соответствуют полученным оценкам в значащих точках, а элементы biT — пропущенным данным. T соответствуют порядковому номеру наблюдения во временном ряду, i — положение значащей точки в окне интерполяции, j — шаг окна интерполяции (совпадает с номером строки), а величина окна выбрана p = 5.
Рисунок 16 — Исходные данные MODIS
Рисунок 17 — Однодневный композит MODIS
Рисунок 18 — Интерполированное изображение MODIS
Рисунок 19 — Исходные данные VIIRS
Рисунок 20 — Однодневный композит VIIRS
Рисунок 21 — Интерполированное изображение VIIRS
На рисунках 16–21 проиллюстрированы ключевые этапы предварительной обработки данных приборов VIIRS и MODIS для отдельно выбранных дней года. Предварительная обработка данных спутниковых наблюдений, описанная выше, позволяет получать изображения, очищенные от влияния шумов, облачного и снежного покрова, теней, а также влияния их краевых эффектов и других мешающих факторов. В результате обработки временных серий и восстановления данных интерполяционным алгоритмом удаётся восстановить информацию, недоступную для прямого использования. При этом получаются восстановленные пространственно однородные изображения и непрерывные временные ряды, что важно при дальнейшем мониторинге и анализе состояния растительного покрова.