4.2.1. Разработка автоматических методов предварительной обработки спутниковых данных среднего и высокого пространственного разрешения с целью фильтрации эффектов влияния мешающих факторов и формирования однородных временных серий данных наблюдений

Для исключения зашумленных и заполнения пропущенных значений во временных сериях дистанционных измерений вегетационных индексов в Отделе технологий спутникового мониторинга ИКИ РАН разработан алгоритм их оконной аппроксимации. Алгоритм предусматривает аппроксимацию временных рядов данных дистанционных измерений на основе полинома второй степени вида с использованием скользящего вдоль оси времени окна, включающего фиксированное количество ненулевых измерений для вычисления коэффициентов полинома. Поиск значений коэффициентов основан на методе наименьших квадратов, подразумевающем минимизацию отклонения значений от аппроксимирующей их кривой. Задача минимизации решается методом деформируемого многогранника (методом Нелдера-Мида), который не использует значения градиентов функции и пригоден для негладких функций, выгодно отличаясь быстродействием и надежностью. Полином второй степени является наиболее подходящим для локальной аппроксимации фенологической динамики растительности, являясь относительно простым и одновременно позволяющим корректно описывать фенологические экстремумы. Вычисление коэффициентов полинома происходит в пределах интервала, занимаемого фиксированным числом валидных измерений, и при этом полином не экстраполируется за пределы окна локализации. Пример скользящей оконной аппроксимации временного ряда значений вегетационного индекса PVI полиномами второй степени с размером окна равным пяти показан на рис.4.2.1.1.
Из рисунка видно, что выбросные значения заметно отличаются от соответствующих значений аппроксимирующей функции, несмотря на то, что они также используются при её построении. Некоторое увеличение размеров окна локализации в случае высокой частоты измерений может приводить к более правильному виду аппроксимирующей кривой и давать больше возможностей для выявления выбросных значений, однако динамика развития растительности на длинных временных интервалах существенно хуже приближается квадратичными функциями. Таким образом, возникает задача поиска оптимального размера окна локализации с учетом особенностей развития растительности в течение вегетационного периода. Поскольку в ряде задач по распознаванию растительности на основе сезонных рядов спутниковых данных используются четырёхдневные интервалы между измерениями с допуском одного потенциально зашумленного значения, а также с учетом скорости роста биомассы, было выбран размер окна локализации в 5 точек, как наиболее оптимальный.

Рис.4.2.1.1. Пример 5-точечной оконной аппроксимации полиномами второй степени. Цветами показаны результаты аппроксимации при различных положениях скользящего окна.

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

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

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

При этом, если   и , то значение текущего элемента заменяется расчетным средним: . Если не выполняется ни одно из этих условий, значение элемента не меняется. Пропущенные значения исходной временной серии также заменяются расчетным средним на основе оценок окружающих локальных полиномов: .

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

Пример полиномиального сглаживания и заполнения пропусков в случае с размером окна локализации равным 5 точкам, количеством проходов равным 3, с использованием временной серии четырехдневных композитных изображений, приведен на рис.4.2.1.2.

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

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

Рис. 4.2.1.2. Результат сглаживания и заполнения пропущенных значений отрезка временного ряда с помощью алгоритма оконной полиномиальной аппроксимации

Рис.4.2.1.3. Снижение зашумленности спутниковых изображений в результате сглаживания с помощью алгоритма оконной полиномиальной аппроксимации, слева – до сглаживания, справа – после сглаживания

В настоящее время достаточно актуальной является задача разработки методов  автоматической обработки спутниковых изображений высокого пространственного разрешения, в том числе направленных на фильтрацию мешающих факторов, таких как облака и тени. К числу наиболее широко используемых спутниковых данных такого рода относятся изображения, получаемые со спутников Landsat-TM/ETM+, с пространственным разрешением около 30 м. Разработанный в отделе технологий спутникового мониторинга ИКИ РАН метод позволяет автоматически детектировать участки облаков и теней на изображениях  Landsat-TM/ETM+.

Метод предусматривает использование предварительно калиброванных данных, выраженных в значениях спектрального коэффициента отражения и радиационной температуры. При этом используются данные измерений в пяти спектральных каналах, включая TM2 (0.60 мкм), TM3 (0.65 мкм), TM4 (0.80 мкм), TM6 (1.65 мкм) и T (11.5 мкм).

Ядром метода детектирования облаков является модернизированный алгоритм ACCA (Automated Cloud Cover Assessment), включающий два этапа. На первом этапе создаётся маска облачности, в которой всем пикселям присваивается флаг облаков и на основе анализа значений каждого пикселя в каналах меняется флаг на “безоблачный” или “сомнительный”. Пиксели, у которых на первом этапе флаг не изменился, относятся к облачным. Для детектирования облаков на первом этапе проводится сравнение коэффициента отражения в канале TM3 и яркостной температуры T с пороговыми  значениями. На втором этапе анализируется распределение температуры облачных пикселей и вычисляются средняя яркостная температура T и среднеквадратичное отклонение σ этого распределения. Все “сомнительные” пиксели, у которых температура выше T+σ, получают флаг “безоблачные”. Оставшиеся “сомнительные” пиксели получают флаг облаков. На рисунке 4.2.1.4 приведен пример детектирования облаков вышеописанным методом.

Рис. 4.2.1.4 Фрагмент изображения Landsat-ETM+ с наличием облакаков (слева) и результат их автоматического детектирования (справа)

Метод  детектирования теней основан на использовании полученной ранее маски облаков и включает в себя два этапа, правый из которых предусматривает построение консервативной, а второй либеральной масок. Консервативная маска может содержать недовывленные участки теней, в то время, как либеральная маска характеризуется некоторой избыточностью. Целесообразность использования того или иного варанта маски теней определеяется  исходя из условий решаемой тематической задачи.

На первом этапе из известных зенитного и азимутального углов Солнца определяется направление, в котором находятся тени и строятся области их нахождения из предположения, что высота облаков не превышает 5 км. В построенных областях предполагаемой тени и вне их строятся гистограммы значений яркости в каналах TM4 и TM5. На основе анализа полученных гистограмм определяются пороги в соответствующих каналах. Всем пикселям со значениями коэффициента отражения  в каналах TM4 и TM5 ниже соотвутствующих порогов присваивается флаг тени.

На втором этапе на основе построенной маски теней анализируется распределение значений расстояния теней от облаков с вычислением значений среднего h и среднеквадратичного отклонения σ.  Затем пиксели, находящиеся на расстояниях от h – 2σ до h+2σ  от облачных пикселей в направление теней, независимо от их яркостных характеристик, получают флаг потенциальной тени. На рисунке 4.2.1.5 показаны примеры создания масок теней после первого и второго этапа обработки соответственно. В первом варианте тени детектируются как после первого этапа, во втором как после второго этапа, который является более надёжным, но приводит к частичной потере данных.

Рис.4.2.1.5. Изображения консервативной (слева) и либеральной (справа) масок теней (белый цвет) от облаков (красный цвет) для врагмента сцены Landsat-5