CC BY
21
Алгоритм выявления случайных искажений в составе сцены на серии разновременных изображений ДЗЗ одной и той же территории
А.М. Белов1, А.Ю. Денисова1 1 Самарский национальный исследовательский университет имени академика С.П. Королёва, 443086, Россия, г. Самара, Московское шоссе, д. 34
Аннотация
Ряд задач обработки разновременных изображений дистанционного зондирования Земли одной и той же территории требует выявления на изображении объектов, не характерных для территории и представляющих собой случайные искажения в составе сцены. К таким искажениям можно отнести облака, тени и другие объекты или результаты воздействия природных явлений, которые перекрывают часть наблюдаемой сцены или существенно меняют регистрируемую яркость объектов в её составе. Случайный характер искажений проявляется в том, что их наличие, расположение, размеры и форма зависят от времени регистрации изображений, т.е. проявляются не на всех снимках из анализируемой серии. В настоящей статье предлагается алгоритм детектирования искажений в составе сцены по серии разновременных изображений дистанционного зондирования Земли. Алгоритм основан на суперпиксельной сегментации изображений и обнаружении аномалий в многомерных потоках данных. Результатом являются маски случайных искажений в составе сцены для каждого из изображений в серии, что позволяет впоследствии в методах комплексирования данных учитывать только релевантные для сцены участки каждого из изображений. Предлагаемый подход отличается универсальностью с точки зрения спектрального и пространственного разрешения анализируемых данных дистанционного зондирования Земли и допускает использование изображений с различной спектральной и пространственной дискретизацией в одном наборе. Качество работы алгоритма оценивалось путём моделирования серии разновременных мультиспектральных изображений с различными параметрами спектральной и пространственной дискретизации при различных условиях облачности и теней от облаков в качестве примера случайных искажений в составе сцены. В результате экспериментальных исследований было показано, что при оптимальном подборе параметров алгоритм обеспечивает точность обнаружения порядка 90 % при ошибке ложного обнаружения порядка 10 %.
Ключевые слова: детектирование случайных искажений в составе сцены, комплексиро-вание изображений дистанционного зондирования Земли, суперпиксельная сегментация изображений, обнаружение аномалий в многомерных потоках данных.
Цитирование: Белов, А.М. Алгоритм выявления случайных искажений в составе сцены на серии разновременных изображений ДЗЗ одной и той же территории / А.М. Белов, А.Ю. Денисова // Компьютерная оптика. - 2019. - Т. 43, № 5. - С. 869-885. - DOI: 10.18287/2412-6179-2019-43-5-869-885.
Введение
В настоящее время одной из основных тенденций в разработке алгоритмов обработки изображений дистанционного зондирования Земли (ДЗЗ), в особенности для анализа экосистем различного уровня, является комплексирование изображений, полученных различными съёмочными системами [1]. Целью ком-плексирования, как правило, является повышение пространственного, спектрального, а также временного разрешения данных ДЗЗ и, следовательно, уменьшение степени неопределённости получаемых с их помощью измерений [2-4]. При этом разновременные изображения ДЗЗ могут содержать случайные искажения в составе сцены, представляющие собой локальные существенные изменения яркости наблюдаемой сцены, не связанные с изменением состояния наблюдаемых объектов, но проистекающие из изменения условий наблюдений территории в момент съёмки. При наличии случайных искажений в составе сцены алгоритмы комплексирования изображений нуждаются в разработке дополнительных инструментов анализа
серии разновременных изображений с целью выявления не характерных для сцены участков изображений (например, облака или тени от облаков) для их устранения из рассмотрения при последующем анализе.
К существующим решениям оценки случайных искажений в составе сцены и их коррекции для серии изображений ДЗЗ можно отнести методы обнаружения облаков и их теней, как самых распространённых примеров таких искажений [5 - 9]. Однако данные методы опираются в основном на использование изображений, полученных одним и тем же сенсором ДЗЗ, и ориентированы на физические свойства спектрального отражения облаков с последующим вычислением теней как побочного продукта, в результате анализа метаданных снимков об угле возвышения солнца над горизонтом в момент съёмки и высоте облачности, получаемой по метеоданным. При этом для обнаружения облаков используется пороговая обработка специальных спектральных индексов [5-7]. Альтернативным подходом является классификация с обучением по данным об облачности, полученным
сенсорами с более низким пространственным разрешением, например, сенсором МОБШ [8].
Определение теней является самостоятельной проблемой особенно при анализе изображений высокого разрешения в пределах городской застройки. Такие алгоритмы работают, как правило, в пределах одного изображения [9, 10]. Выделяют два типа детекторов теней на основе анализа яркостей пикселей и на основе моделирования теней с учётом формы и расположения объектов сцены [11]. Первый подход использует предположение о том, что затенённые области имеют более низкую яркость, чем аналогичные не затенённые области, например, в [12] для реализации такого подхода используется медианная фильтрация ряда разновременных значений яркости, получаемых для каждого пикселя. Второй подход требует приближённого определения геометрии объектов сцены, тогда с использованием параметров съёмки (угол возвышения солнца и т.п.) может быть рассчитано положение теней на изображении. В случае с тенями от облаков второй подход также требует задания ряда метеорологических параметров.
Приведённый обзор показывает, что существующие алгоритмы характеризуются узкой спецификой выделяемых случайных искажений в составе сцены и ориентированы на определённые спектральные диапазоны входных данных ДЗЗ, что лишает их универсальности при анализе данных, полученных различными сенсорами ДЗЗ.
В рамках настоящей статьи предлагается более универсальное решение, позволяющее оценить случайные искажения в составе сцены в более широком смысле, нежели облачность и тени, а также более универсальное в отношении характеристик спектрального и пространственного разрешения анализируемых данных. Кроме того, предлагаемый алгоритм не требует привлечения дополнительных априорных данных в виде метаданных снимков или данных метеорологических наблюдений.
Основу предлагаемого алгоритма составляют методы суперпиксельной сегментации изображений, кластеризации и обнаружения аномалий в потоках многомерных неструктурированных данных. При этом детектируемые искажения в составе сцены интерпретируются как временные и пространственно-спектральные аномалии в составе сцены. Алгоритм анализирует стек разновременных изображений, приведённых к единой сетке спектрально-пространственных координат. Суперпиксельная сегментация используется для выделения пространственно локализованных групп пикселей со схожим спектрально-временным профилем на всех изображениях. В рамках каждой такой группы пикселей выделяются спектрально-временные кластеры без учёта пространственных связей. И, наконец, для полученных спектрально-временных кластеров оценивается мера аномальности отдельных спектральных фрагментов, соответствующих различным моментам времени. Результирующие спектральные аномалии в рамках
спектрально-временных кластеров отображаются в виде маски искажений в составе сцены для соответствующих изображений в серии.
Для оценки меры аномальности в рекомендуемой реализации используется алгоритм Local Outlier Factor (LOF) [13]. Мотивация применения алгоритма LOF основана на результатах сравнительного исследования [14], которое показало, что LOF демонстрирует высокое качество обнаружения аномалий на большинстве источников данных, в отличие от других классических методов машинного обучения.
Для проверки качества обнаружения случайных искажений в составе сцены было произведено экспериментальное исследование на серии модельных разновременных изображений двух мультиспектральных сенсоров с различными характеристиками пространственной и спектральной дискретизации. В качестве случайных искажений в составе сцены использовались облака и тени от облаков, искусственно встроенные в ряд снимков. В результате были выбраны рекомендуемые параметры алгоритма и исследовано качество его работы при различных характеристиках входного набора данных. Результаты экспериментов показали, что достижимая при оптимальных параметрах алгоритма средняя точность обнаружения искажений в составе сцены составляет порядка 87-90 %.
1. Постановка задачи
Будем полагать, что идеальное дискретное представление наблюдаемой сцены X содержит L спектральных каналов X = (Xi, ..., Xx, ..., XL), определенных для диапазонов длин волн [Çi(l), Ç2(X)], объединение которых соответствует полному наблюдаемому
диапазону спектра U^[Çi (X),Ç2(X)]=[wmm,"max].
Без ограничения общности положим, что изображение X определено для дискретной сетки координат m1, m2 с шагом дискретизации T = 1, которая определяет эталонное расположение и размер кадра в виде прямоугольного фрагмента 0 < m1 < M1, 0 < m2 <M2.
Пусть сцена наблюдается с помощью K различных систем ДЗЗ в моменты времени tkÇ6[tmin, tmax], k = 1, ..., K, Ç = 1, ..., Çk, где Çk - количество кадров, полученных k-й системой съёмки. Изображение, сформированное k-й системой, будем обозначать Y*ç и будем полагать, что оно имеет шаг пространственной дискретизации Tk > T и содержит Lk спектральных каналов Ykç, 1 < l < Lu.
Состав сцены, а значит, и соответствующее ей идеальное дискретное изображение X будем полагать неизменным в течение всего периода наблюдений
[tmin, tmax].
При отсутствии случайных искажений в составе сцены будем описывать связь между наблюдаемыми изображениями Ykç и идеальным дискретным изображением сцены X в виде:
Ykçi = DkHiFkçXL=1 ^ы (X)XX + Vkçi, k = 1,...,K, l = 1,...,Lk, Ç = 1,...,Çk.
(1)
Модель наблюдения (1) объединяет модель спектральной дискретизации, основанную на модели описанной в [15], и стандартную модель формирования изображений оптическими системами в каждом канале при многокадровой съёмке, описанную в [16]. Модель (1) определяет спектральные и пространственные искажения регистрируемых изображений сцены следующим образом.
Коэффициенты wkl (X) задают пересчёт значений спектральной сигнатуры каждого пикселя (ш\, т2) изображения X в значения наблюдаемого к-й системой 1-го спектрального канала. Каждый канал 1 определяется своей непрерывной функцией спектрального отклика Жк1 (и). Для формирования единой модели для всех К систем принимается предположение, что идеальному дискретному спектральному описанию сцены в виде набора каналов (ХЛ }£=1 соответствует непрерывное кусочно-постоянное спектральное представление, такое, что на каждом интервале [^1(1), <^2(Х)] значение непрерывного спектрального отклика полагается равным Хх. В этом случае коэффициенты пересчёта wкl (X) согласно традиционной модели спектральной дискретизации из [15] можно определить в виде:
(X) = ^ (и) £ (и)Аи . (2)
Операторы Дк%, Нк и Вк определяют соответственно геометрическое преобразование эталонного кадра в наблюдаемый, размытие, вызываемое дефокусиро-вой оптической системы, и пространственную дискретизацию соответственно. В отличие от модели, рассмотренной в [16], в настоящей статье описывает разницу в геометрии кадра не в связи с движением камеры, а в результате остаточного смещения кадров наблюдаемых изображений из-за погрешности геопривязки и в связи с различным размером кадра из-за отличий в шагах дискретизации наблюдаемых изображений Тк. Оператор Вк определяется как равномерная дискретизация с усреднением с шагом дискретизации Тк, что даёт возможность учесть размытие, вызываемое усреднением по апертуре детекторного элемента изображающей системы. Величина Ук%1 задаёт статистически независимый от сигнала белый шум. Все операторы, входящие в модель (1), полагаются удовлетворяющими свойствам линейности и ограниченности.
Для удобства описания обозначим далее оператор спектральной дискретизации к-й системы съёмки за Лк. Очевидно, что Лк является линейным и ограниченным оператором как линейная комбинация спектральных каналов изображения X. Тогда модель наблюдения (1) примет вид:
ВкНкДКЛ кХ + ^, (3)
1 = 1,...,Ь, к = 1....,К. ()
При наличии случайных искажений в составе сцены будем полагать, что на части изображения X в момент времени существенно меняется яркость
пикселей в результате перекрытия непрозрачными объектами (облака) или в результате некоторого явления (тень от другого объекта или облака). Пространственное положение искажений в составе сцены будем описывать бинарной маской искажений Бка: Бка(ш\,ш2) = 1, если в точке ш\,ш2 есть искажения в составе сцены, и 5^(ш1, ш2) = 0 иначе.
Будем называть изображением искажений в составе сцены для момента времени такое изображение Ск%, что в точке сцены ш\, ш2 с искажениями в составе сцены Скфп\, ш2) ФX(ш\, шг), иначе Скфп\, ш2) равно нулевому вектору.
Тогда модель наблюдения разновременных изображений при наличии случайных искажений в составе сцены будет иметь вид: = ВкНкДк, Л ^' + V*,
X ' = [(1 - ) + БКСК ], (4)
, = 1,...,, к = 1,...,К.
Задача обнаружения случайных искажений в составе сцены заключается в оценке масок для каждого из наблюдаемых изображений Ук% при неизвестном идеальном изображении сцены X.
Связь с задачей обнаружения аномалий
Так как изображение сцены без учёта случайных искажений в составе сцены X полагается неизменным или меняющимся незначительно в течение всего периода наблюдений [4яя, ш*], то неизвестное исходное изображение можно рассматривать как сложный фон для объектов, изменяющих сцену, а сигнатуры этих объектов - как аномальные сигнатуры во временном ряду каждого пикселя изображения.
Переведём наблюдаемые изображения Ук% в единую систему спектрально-пространственных координат идеального кадра. Для этого выполним преобразование изображений Ук% в требуемое представление 2к% с помощью пространственной интерполяции 1к с шагом Т/Тк, пространственного преобразования к положению и форме эталонного кадра Дк-1, а также оператора Qk, определяющего преобразование Ьк наблюдаемых каналов в Ь каналов, соответствующих каналам сцены X:
^ка, = QkFk-1IkYk,. (5)
Операторы 1к, и Qk задаются как линейные и
ограниченные для согласования с моделью (1). Тогда, имея в виду (4) и (5), можно записать:
^ =(1 - Бка )А^ + + V ^ , (6)
где Ак, = QkF-lIkВkHkFщ Л kX, В к, = QkFk-1IkВkHkFki;ЛкСка, V = QkFk-IkVkа.
Таким образом, для последовательности разновременных изображений в каждой точке с пространственными координатами составляющая Кк% определяет ожидаемое истинное значение сигнатуры пикселя, а Вка - аномальное. Анализируя временные ряды для каждой точки (ш\, ш2), можно получить значения маски искажений в составе сцены как:
Ёк§(ш1, т2) = 1, если 2к§(ш1, т2) является аномалией во временном ряду
{¿к§(т 1,Ш2))к=1,...л , §=1,...,§ к
и (ш1, ш2) = 0 иначе.
Очевидно, что пространственные и спектральные искажения, связанные с регистрацией изображений Ук§ и их преобразованием в представление 1к%, должны быть существенно меньше по норме, чем разность образов пикселей для участков с искажениями в составе сцены и без искажений в составе сцены:
В^, (ш!, Ш2) - А^, (ш1, Ш2 )|| >
>||Ак,1;1 (ш1, ш2 )-Ак2?2 (шь ш2 )||, (7)
к * к2,§ 1 * §2.
Условие (7) определяет ограничение на возможную область решения поставленной задачи, а именно: детектируемые искажения в составе сцены должны иметь спектральную сигнатуру, существенно отличную от истинного значения спектральной сигнатуры в анализируемой точке (ш1, ш2) и её окрестности. Кроме того, очевидно, что преобразование QkFk-§1Ik
следует выбирать таким образом, чтобы оно как можно меньше увеличивало мощность шума, так как очевидно, что детектирование изменений, сопоставимых по амплитуде и мощности с шумовой составляющей, бессмысленно.
При соблюдении обозначенных выше условий, являющихся естественными ограничениями, определяемыми свойствами модели наблюдения, задача обнаружения случайных искажений в составе сцены может быть сведена к задаче поиска аномалий в последовательности разновременных наблюдений 1к% в каждой пространственной точке (ш1, ш2). Однако, так как обеспечить устойчивый результат для одной точки изображения на коротких временных сериях достаточно сложно, то имеет смысл перейти к решению задачи поиска аномалий для пространственно и спектрально однородной группы пикселей.
Предположим, что для набора разновременных изображений
^к§ }к =1,..., К §=1,..., §к
известно пространственное разбиение множества [0, М1 - 1]х[0, М2 - 1] координат пикселей на спектрально и пространственно однородные подмножества {Оу}у=1!...!г:
Оу = {(шьш2): У(шьш2), (ш 1,ш'2) е л
л(к1 * к2 V * §2 )
^ р (щ (ш1, ш2), 2к2§2 (ш 1, ш '2)) +
+^р((/и1,ш2),(ш'1,ш'2)) <г|,
где р - евклидово расстояние, ^ - весовой коэффициент, г - некоторая константа. Обработка пикселей всех разновременных изображений, принадлежащих
одному подмножеству Оу, вместо попиксельной обработки, позволяет увеличить количество отсчётов, используемых для решения задачи поиска аномалий, а также учесть такое свойство изображений, как пространственная коррелированность.
Окончательно получаем, что решаемая в настоящей статье задача отыскания искажений в составе сцены сводится к поиску аномалий во временных рядах пространственно коррелированных пикселей
(шь ш2 )=1,...,к
((ш ,ш2
и в отображении найденных аномалий на множество бинарных изображений искажений в составе сцены §к§, к = 1, ..., К, § = 1, ..., §к. При этом полагается выполненным условие (7). Изображение идеального дискретного представления сцены X полагается неизвестным. Операторы Qk, Дк-1,1к считаются заданными и согласованными с моделью наблюдения изображений Ук§. Оценка модели наблюдений и оптимальный выбор операторов интерполяции Qk, 1к выходят за рамки настоящего исследования.
2. Предлагаемый алгоритм
Общее описание
Схема предлагаемого алгоритма приведена на рис. 1. Исходными данными алгоритма является набор разновременных изображений одной и той же территории Ук§, к = 1, ..., К, § = 1, ..., §к, сформированных К различными системами ДЗЗ и совмещённых с точностью до погрешности геопривязки. Далее, для обозначения общего количества изображений во входном наборе данных будем использовать обозначение: Н = ХК=1§к .
Модель наблюдения полагается известной и определённой относительно параметров эталонного кадра наблюдаемой сцены, т.е. шага пространственной дискретизации Т и характеристик для Ь спектральных каналов.
Выходными данными является набор бинарных изображений оценок искажений в составе сцены §], ] = 1, ..., Н в системе пространственных координат эталонного кадра. Заметим, что получаемые маски имеют большее пространственное разрешение, чем исходные сопоставляемые снимки, поскольку Т < Тк, к = 1, ..., К.
В общем случае предлагаемый алгоритм поиска случайных искажений в составе сцены может быть сформулирован в следующем виде:
1. Исходные изображения Ук§ Ук§, к = 1, ..., К, = 1, ... , к приводятся к системе пространственно-
спектральных координат эталонного кадра согласно (5). В результате формируется набор изображений 1к§, Ук§, к = 1, ..., К, § = 1, ..., §к, согласованных по размеру и количеству спектральных каналов.
2. Для стека изображений 2к§, имеющего размер М1 х М2 х НЬ, выполняется суперпиксельная сегментация. В результате получается единое для всех изоб-
ражений разбиение на спектрально и пространственно однородные подмножества координат пикселей Оу}у=1,...,г. Далее будем называть Оу суперпикселями.
3. Для каждого суперпикселя независимо решается задача поиска аномалий:
3.1. Для этого последовательность пикселей
{гк (ш1, ш2 )}^==1:::кк ,
(ш ,ш2 )еОу
принадлежащих О- на всех изображениях 1к%, преобразуется в последовательность спектрально-временных пикселей {г/ (шь ш2)}(ш1> ш2)ЕО,. Каждый вектор г/ (ш1, ш2) е ЯЕЬх1 является конкатенацией Ь-мерных
векторов-пикселей гука (ш\, ш2), к = 1, ..., К, , = 1, ..., в точке (шь ш2).
3.2. Для учёта спектрально-временной корреляции пикселей выполняется кластеризация последовательности {г/ (шь шг)}^,,^)^ на О спектрально-временных кластеров. Количество кластеров О определяется по формуле:
О = Е/ Е, (9)
где Е - параметр алгоритма. Полученные центры спектрально-временных кластеров обозначим как
{г: (у)и.., О.
к=1
Стек разновременных изображений от К систем ДЗЗ
I I I I I I I I I I |Д
Последовательность Ь-мерных центров спектральных кластеров для К снимков в суперпикселе Су
Приведение к единым
спектральным и пространственным параметрам кадра
Суперпиксельная сегментация
Получение центров спектральных кластеров для каждого снимка из спектрально-временных центров кластеров
Спектрально-
временная кластеризация пикселей из С
Вычисление меры аномальности
центров спектральных кластеров в суперпикселе ву
Отображение аномальных
суперпикселей или их аномальных подкластеров на результирующие маски искажений в составе сцены
Разбиение на пространственно однородные подмножества Gy.gr!.....Г
Бщ, к=1,...,К, х=1....,\к
3.3. Для перехода к оценке аномалий последовательность спектрально-временных центров кластеров {г: (/)}„=! ..., О преобразуется в последовательность спектральных центров кластеров
{М}о=1,.О
0=1,...=
Рис. 1. Схема алгоритма
{Фо/(/)}»=!,...,о .
без учёта временной составляющей, где / - общий индекс изображения гк в наборе, определяемый для фиксированного к по формуле:
/=Уа к+а, а=1,...,а к.
К=1
(10)
Для этого каждый вектор г» (у) е ЯЕЬх1 разбивается на Е векторов г'»;(у) е ЯЬх1, соответствующих отдельным изображениям гк . 3.4. Для последовательности
{г»; (у)}»=1,..., О
выполняется оценка мер аномальности
Полученные значения фо(у) используются для принятия решения об отнесении суперпикселя О- частично или полностью к искажениям в составе сцены (аномалиям) на каждом из Е анализируемых изображений.
3.5. Для начала проверяется статистическая значимость отличия в среднем подпоследовательности мер аномальности ф»(у), о = 1, ..., О для изображения / от всей последовательности
{Фо;(у)}»=1,...,О .
Проверка производится путём использования статистического теста Стьюдента для двух выборок с критерием значимости 0,05 [17].
Если отличие в среднем для данных последовательностей значимо, то для изображения с индексом / суперпиксель Оу целиком относится к области искажений в составе сцены, т.е. (шьш2) = 1, (ш\,ш2) е Оу, где индексы к, связаны с индексом/ соотношением (10).
3.6. В случае, если суперпиксель Оу для изображения ] частично затронут искажениями в составе сцены, тест, выполняемый на шаге 3.1, вероятнее всего не пройдёт. Для таких случаев необходимо определить, какие из центров спектральных кластеров суперпикселя Оу являются аномальными относительно всей рассматриваемой последовательности
{^ (у)ь=1,.., о,
т.е. содержат изменения в составе сцены. Для этого меры аномальности {Ф» (У)}»=1,..., о
подвергаются пороговой обработке:
(у) =
[1, ф„ (у)>ф(у), [о, Фо(у)<ф(у),
(11)
где ij) (у) - у-квантиль распределения значений мер аномальности в последовательности
(Фо/ (У)}о=1,...,0 •
Если для изображения с индексом j выполнено:
£ и„ (у)>ю0 , (12)
0=1
где гае [0,1] - параметр алгоритма, то пересечение суперпикселя Gy с областью искажений в составе сцены полагается существенным и пиксели (m1, m2), принадлежащие кластерам е, маркируются как содержащие искажения в составе сцены: §к% (m1, m2) = 1 (связь индексов Щ и j задаётся соотношением (10)). Если (12) не выполнено, полагается, что существенных искажений в составе сцены для суперпикселя Gy на изображении j не обнаружено.
4. Если на шаге 3.5 и 3.6 ни одно из условий аномальности для суперпикселя Gy на изображении j не было выполнено, то значение Sk% (m1, m2) = 0 для всех (m1, m2)eGy (связь индексов Щ и j задаётся соотношением (10)).
Рекомендуемая реализация алгоритма
Изложенный выше алгоритм не конкретизирует методы и алгоритмы суперпиксельной сегментации, кластеризации и вычисления мер аномальности, таким образом, подразумевая возможные вариации реализации данных этапов. Тем не менее, основываясь на наблюдениях авторов настоящей статьи и результатах предварительных экспериментов, можно выделить следующие предпочтительные методы и алгоритмы.
Суперпиксельную сегментацию на шаге 3.1 предлагается производить алгоритмом Simple Linear Iterative Clustering (SLIC), предложенным в [18]. Метод SLIC предназначен для выделения на изображении пространственно и спектрально коррелированных областей - суперпикселей, и, согласно обзору [19], обладает минимальным временем работы и ошибкой пересегментации в классе аналогичных методов. Основу SLIC составляет алгоритм кластеризации k-means [20] с мерой близости, аналогичной расстоянию (8). В ори-
гинале SLIC использует переход в цветовое пространство LAB [21] для расчёта спектральной составляющей меры (8). Однако для изображений ДЗЗ переход к пространству LAB может быть в общем случае не определён, так как набор спектральных каналов может существенно отличаться от RGB-представления. Поэтому в настоящей работе при реализации SLIC используется полное спектральное евклидово расстояние между пикселями. Весовой коэффициент ^ в выражении (8) является параметром алгоритма и определяется из соображений формирования пространственно компактных кластеров.
Для выделения спектрально-временных кластеров на шаге 3.2 предлагается использовать алгоритм k-means [20], обладающий высокой скоростью работы и простой реализацией.
Для расчёта меры аномальности рекомендуется использовать метод Local Outlier Factor (LOF) [13], который согласно сравнительному исследованию в [14] обеспечивает высокое качество обнаружения аномалий на большинстве наборов данных. Мера LOF показывает, насколько изолированно расположен в пространстве признаков проверяемый вектор признаков. Метод LOF имеет один основной параметр - количество ближайших соседей P. Если тестируемый вектор признаков в среднем расположен существенно дальше от своих P ближайших соседей, чем каждый из них от своих P ближайших соседей, то данный вектор признаков рассматривается как аномалия. По сути, LOF оценивает аномалии на основе свойств локальной плотности распределения векторов признаков. Подробности метода LOF изложены в работе [13].
Реализация предложенного алгоритма с помощью рекомендуемых методов требует учёта следующих нюансов.
Во-первых, метод SLIC не ограничивает минимальное количество точек в суперпикселе. Поэтому при малом размере суперпикселя:
\ву\< 3O (13)
требуется адаптивно снижать количество кластеров O:
O = ЦСу|/3] . (14)
Выражение (14) определяет условие наличия по крайней мере трёх точек в кластере в случае, если размер суперпикселя недостаточен для формирования количества кластеров O, определяемого формулой (9).
Во-вторых, значение параметра P для расчёта меры LOF зависит от мощности последовательности {Фо (у)}»=1,..., о .
Чем меньше длина последовательности, тем, очевидно, меньше должно быть P, так как в противном случае теряется свойство локальности отличий в плотности распределения. Поэтому предлагается использовать два значения параметра P: P1 в общем случае и P2 в случае, когда выполняется (13). И, наконец, если размер суперпикселя настолько мал, что P > OS, то можно положить:
P = |OH/ 3J .
(15)
Выражения (13) - (15) были определены в результате экспериментальных исследований и позволяют контролируемо перенастраивать метод LOF в зависимости от размера анализируемого суперпикселя.
Окончательно параметры рекомендуемой реализации предложенного алгоритма имеют вид:
1) весовой коэффициент учёта пространственной составляющей расстояния для суперпиксельной сегментации SLIC - ^е [0, 100],
2) количество суперпикселей Г, выделяемых на изображении, определяет также средний размер суперпикселя как M\ M2 / Г,
3) параметры метода LOF для крупных и малых суперпикселей - P1 и P2 соответственно,
4) параметр уе[0,1] для определения порога аномальности при частичном пересечении суперпикселя с областью искажений в составе сцены,
5) длина E последовательности, анализируемой на наличие аномалий методом LOF для крупных суперпикселей, используется для вычисления количества кластеров в формуле (9),
6) юе[0,1] - допустимая доля пропускаемых аномальных кластеров при частичном пересечении суперпикселя с областью искажений в составе сцены.
3. Схема экспериментальных исследований
Показатели качества оценки искажений в составе сцены
Исследование предложенного алгоритма производилось путём моделирования серии разновременных мультиспектральных изображений, полученных различными системами ДЗЗ с искажениями в составе сцены и без них. Качество работы алгоритма оценивалось путём сравнения эталонных масок искажений в составе сцены
{SkÇ }k=1,..., K Ç=1,..., Çk
с их оценками
{SkÇ }k=1,..., K .
Ç=1,..., Çk
Для изображений с искажениями в составе сцены оценивались средняя ошибка ложного обнаружения p1 и средняя ошибка пропуска искажений p2. Для изображений без искажений оценивалась только средняя ошибка ложного обнаружения p'1.
Выражения, используемые для расчёта ошибок, приведены ниже. Для удобства записи полагается, что изображения с искажениями в составе сцены имеют индексы j = 1, ..., S, а изображения без искажений - j = Si+1, ..., H - Si соответственно:
H, M1-1M2-1 . .
( m2 ) - Sj (m, m ) j Sj (mb m2 )
P1 =
j=1 mt\ =0 m2 =0
M1M2S1
H1 M1-1M2-1
SSS(Sj (m1, m2 ) - S j (m1, m2 )) S j (mb m2 )
_ j=1 m =0 m2=0 p2 =-
H-Hi Mi-1 M2-1
S ££(( (1,m2)-Sj (m1,m2))) (m1,m2)
M1M2H1
,(16)
I _ j=H1 m1=0 m2=0
p1 -
M1M2 (H-H1 )
где Н1 - количество изображений с искажениями в составе сцены.
Моделирование изображений с искажениями в составе сцены
В качестве случайных искажений в составе сцены в экспериментах рассматривались облака и их тени, поскольку они встречаются во многих практических случаях.
Эталонное X и наблюдаемые Ук%, к = 1, ..., К, % = 1, ..., ^к изображения сцены формировались для одного и того же безоблачного участка гиперспектрального снимка Хн, полученного с помощью сенсора ЛУТОК [22]. При этом для идеального дискретного представления Х использовалась только спектральная дискретизация, а для наблюдаемых изображений при отсутствии в них изменений использовалась модель наблюдения (1) с гиперспектральным изображением Хн в качестве входного изображения.
Для моделирования изображений с искажениями в составе сцены перед применением модели наблюдений
(1) в гиперспектральное изображение Хн встраивались облака и их тени. Обозначим гиперспектральное изображение с искажениями в составе сцены Х'нк%, где индексы к показывают, что данное изображение будет использовано в качестве исходного в модели (1) для получения наблюдаемого изображения У^.
Используемый в настоящей статье алгоритм имитации облачности подходит для моделирования непрозрачных кучевых облаков, а принцип формирования теней основан на методе, описанном в [23].
При моделировании предполагалось, что гиперспектральная сигнатура облака постоянна в пределах облака и равна с*. Пространственное расположение облака определялось маской облачных пикселей БС1'1. Для формирования теней для каждого изображения Х'нк% задавались угол возвышения солнца над горизонтом аь (к, | ), азимутальный угол аА (к,| ) и высота нижней границы облака Ис.
Используемый алгоритм формирования изображений Х'нк% имеет вид:
1. Вычислить угол направления на солнце в зените:
а2 (к) = (90-аь (к))/180. (17)
2. Сформировать маски теней от облаков Б^"1™ как результат пространственного сдвига маски облаков с масштабированием:
Shadow (, ) = (m1 + Дь m1 + Д2), Ai = tan(az (k,l)) hc • cos( (k,l)), Д2 = tan (a z (k, £))• hc • sin ((k, l)),
(18)
где (ш'1,ш'2) - координаты тени от облачного пикселя (ш1, ш2), кс - высота нижней границы облака.
3. Сформировать изображение с облаками и тенями:
X V = (( - VcsShdow) (1 - Sit") -
(19)
где ц = 0,25 - коэффициент, задающий уровень снижения яркости пикселей исходного изображения тенью.
4. Сформировать результирующую маску искажений в составе сцены &§ как объединение множеств облачных пикселей и пикселей, попавших в тень от облака:
О _ Oshadoy
| S cloud
(20)
где | - поэлементное логическое «или».
Сигнатура облака cs формировалась как средняя сигнатура участка кучевого облака с того же изображения AVIRIS, что и для формирования эталонного и наблюдаемого изображений. График использованной в эксперименте сигнатуры облака приведён на рис. 2.
Для формирования реалистичных границ облаков Schud применялась пороговая обработка выборочных фрагментов изображений 9 канала Landsat 8 (Cirrus band). Согласно [24], в этом канале (1375 нм) высокое поглощение парами воды, в результате фон под облаками является тёмным и может быть удалён путём пороговой обработки.
12000
500
700
900 1100 1300 Длина волны, нм
Рис. 2. Спектральная сигнатура облака, использованная в эксперименте
Параметры пороговой обработки для каждого изображения Х'нк§ подбирались так, чтобы обеспечить заданное количество пикселей с искажениями в составе сцены (т.е. с учётом теней). Параметры, определяющие положение теней от облаков, задавались случайным образом в пределах следующих интервалов: аке [60, 80], а^е [15, 350], ^е [3100, 5900].
Примеры изображения сцены без искажений Хн и изображения сцены с искажениями Х'нк§ приведены на рис. 3.
Параметры модели наблюдения
В экспериментах моделировались изображения для двух систем ДЗЗ с различными параметрами про-
странственной и спектральной дискретизации. Общие характеристики эталонного изображения сцены и наблюдаемых изображений приведены в табл. 1.
Рис. 3. Пример изображения без случайных искажений в составе сцены (а), с искажениями в виде облака и тени от облака (б)
Табл. 1. Общие параметры эталонного и наблюдаемых изображений
Параметр Эталон Система 1 Система 2
Размер 256x256 128x128 64x64
Шаг простран-
ственной дискрети- 1 2 4
зации
Количество каналов 16 6 4
Спектральная дискретизация эталонного и наблюдаемых изображений сцены производилась с помощью гауссовых функций спектрального отклика:
W, (и) = 571 (2я) exp {-0,5 (и - u,°) 2/s?}, (21)
где и0 - центральная длина волны, а ширина спектральных каналов (FWHM [25]) определяется как . FWHM используется в документации систем ДЗЗ для описания спектрального состава снимков и позволяет произвести моделирование в соответствии с параметрами реальных мультиспектральных систем ДЗЗ.
В табл. 2 представлены параметры и0 и FWHM для эталонного и наблюдаемых изображений. Система 1 имеет параметры спектральных каналов сенсора Геотон (Ресурс-П) [26], а система 2 имеет параметры сенсора SPOT-7 [27].
Пространственные искажения наблюдаемых изображений определялись в соответствии с моделью (1) и были заданы следующим образом:
X^ (1, m2, /) = Fk% (( (1, m2, /)) =
= (1 + Xiy, m2 + Xk;2, l),
(mm,m2,/) = Hk ( ((,m2,/)) =
3« (23)
= X V(T1, T2) Xi?(w1 -Tb Ш2 -T2, l),
TbT2 =-3nt
v(x1, T2) = A exp (-0,5 -ст-2 (т2 + т2) ) , (24)
XD (1,Ш2,l) = Dk1 ( (mbm.2,l)) =
^/2] [n/2-] (25)
= X X xh (+T1,m2+T1,1),
T1 =-\Ttj 21 T2=-\Ttl 21
(22)
(,«2,l) = Dk2 ( (m2,l)) =
= ( m2 , l )) =„Л , (26)
m=n2rt
где Ха; = Хн для изображений без искажений в составе сцены и Ха; = ХНа; иначе, Fa; - сдвиг кадра на ха; 1, 2 пикселей, На - Гауссово размытие, Dai - усреднение по апертуре детекторного элемента, Dk2 - децимация. Операторы Dai и Dk2 задают дискретизацию с усреднением по площади пикселя наблюдаемого изображения.
Табл. 2. Параметры спектральной дискретизации эталонного и наблюдаемых изображений
Номер канала Эталон Система 1 Система 2
Ul FWHM Ul FWHM Ul FWHM
1 462,75 28,82 485 70 490 70
2 491,90 28,75 560 80 560 60
3 521,08 28,71 645 70 660 70
4 550,30 28,68 685 30 825 130
5 579,55 28,68 715 30 - -
6 608,85 28,70 750 100 - -
7 638,19 28,73 - - - -
8 665,07 29,46 - - - -
9 694,47 29,31 - - - -
10 723,83 29,25 - - - -
11 753,14 29,27 - - - -
12 782,41 29,38 - - - -
13 811,62 29,58 - - - -
14 840,79 29,87 - - - -
15 869,92 30,24 - - - -
16 898,99 30,71 - - - -
Параметры i, ^ 2 задавались случайным образом в пределах интервалов [-1, 1] и [-2, 2] для системы 1 и системы 2 соответственно. Импульсная характеристика (24) определяет Гауссово размытие, причем константа A обеспечивает выполнение условия нормировки, а радиус пятна размытия для системы 1 и системы 2 был равен ст1 = 2, ст2 = 4 соответственно.
Количество изображений Ye, для каждой из систем было выбрано: Е 1 = 4 и Е 2 = 16, - исходя из условия оптимальности восстановления изображений по множеству кадров с повышением разрешения [28]: Е k = T2. Ориентация на алгоритмы повышения разрешения изображений связана с дальнейшими планами по применению предложенного алгоритма для маскирования искажений в составе сцены в процессе восстановления изображения с повышением разрешения.
Программная реализация
Программная реализация предложенного алгоритма и его экспериментальное исследование были выполнены в среде MATLAB 2017 [29]. Алгоритм LOF использовался в виде реализации, описанной в статье [30]. Для реализации алгоритма SLIC была модифицирована версия алгоритма, представленная в статье [19], суть модификации изложена в подразделе «Рекомендуемая реализация алгоритма».
4. Результаты экспериментальных исследований
Выбор оптимальных параметров алгоритма
Настройка предложенного алгоритма была произведена на базовом наборе из 20 тестовых изображений, содержащих 50 % изображений с искажениями в составе сцены. Доля пикселей с искажениями в составе сцены была равна 10 % от общего количества пикселей на изображении.
Рекомендуемые значения параметров алгоритма приведены в табл. 3. Результаты для других значений параметров не приводятся ввиду большого количества возможных вариантов, однако выбранные значения обеспечивают оптимальное соотношение между ошибками ложного обнаружения и пропуска искажений при сохранении обеих ошибок на уровне 10 %.
Табл. 3. Рекомендуемые параметры алгоритма
Параметр Тестовый диапазон Рекомендуемые значения
Л от 10 до 100 с шагом 10 60
Г от 500 до 3000 с шагом 500 2000
P1 от 10 до 40 с шагом 5 20
P2 от 10 до 20 с шагом 5 10
V от 0,05 до 0,3 с шагом 0,05 0,1
E 32, 64, 128, 256 64
ю 0, 0,25, 0,5 0
Пример результатов, полученных при рекомендуемых параметрах, приведён на рис. 4. Для данного случая на изображениях с искажениями в составе сцены средняя ошибка ложного обнаружения составила pi = 0,08, ошибка пропуска искажений - pi = 0,10. Для изображений без искажений в составе сцены средняя ошибка ложного обнаружения была равна
p'1 = 0,01.
Наибольшее влияние на скорость работы алгоритма оказывает параметр E - длина последовательности, используемой для анализа аномалий алгоритмом LOF. Рост трудоёмкости алгоритма LOF с увеличением E связан с использованием поиска ближайших соседей для оценки меры аномальности. Рекомендуемое значение E было выбрано из соображений удобства проведения экспериментов и обеспечивало время обработки всего набора данных менее, чем за 2 минуты.
Параметры Pi и Pi, определяющие количество анализируемых ближайших соседей при расчёте меры аномальности, также зависят от E и были подобраны экспериментально для обеспечения значений ошибок пропуска искажений и ложного обнаружения на уровне 10 %.
Параметр у, используемый для определения порога аномальности при частичном пересечении суперпикселя с областью искажений, оказывает наиболее существенное влияние на ошибки pi и pi среди всех параметров алгоритма. Увеличение у приводит к сокращению ошибки пропуска искажений и росту ошибки ложного обнаружения.
Рис. 4. Пример результатов работы алгоритма при рекомендуемых параметрах (исходные изображения приведены к одному размеру для удобства восприятия): для системы 1: а) исходное изображение, б) эталонная маска искажений в составе сцены, в) результат оценки искажений; для системы 2: г) исходное изображение, д) эталонная маска искажений в составе сцены, е) результат оценки искажений
Весовой коэффициент ^ в выражении (8) был определён как доставляющий наилучшее визуальное качество для суперпиксельной сегментации в смысле пространственной и спектральной однородности получаемых областей.
Количество суперпикселей Г определяет средний размер суперпикселя. Очевидно, что чем больше суперпиксель, тем больше статистики учитывается при оценке меры аномальности и ошибка пропуска искажений р>2 сокращается, но возможно увеличение ошибки ложного обнаружения р\, и наоборот. Параметр ю имеет смысл только для очень крупных суперпикселей (при Г = 500), и в большинстве случаев может быть положен равным нулю.
Проведённые эксперименты показали, что рекомендуемые значения для Е, Р1, Р2, ^ и ю можно считать фиксированными. Поэтому дальнейшее исследование зависимости качества работы алгоритма от объёма анализируемых данных и соотношения количества изображений с искажениями в составе сцены и без искажений проводилось для параметров Г и у.
Зависимость качества работы алгоритма от общего количества изображений в наборе данных
Одним из важнейших факторов, влияющих на качество работы алгоритма, является общий объём анализируемых данных. Поэтому был проведён эксперимент по оценке зависимости качества работы алгоритма от общего количества изображений в наборе.
Были рассмотрены выборки из исходного набора данных, содержащие от 4 до 20 изображений для обеих рассматриваемых систем. Количество изображений с искажениями и без было 1:1. В качестве изменяемого параметра алгоритма выступало количество выделяемых суперпикселей Г. Для остальных параметров были использованы рекомендуемые значения, указанные в табл. 3. Результаты эксперимента приведены на рис. 5.
С уменьшением количества изображений до 6 ошибка ложного обнаружения p\ растёт постепенно с 9 до \4 % при любых значениях параметра Г. При этом ошибка пропуска искажений p2 остается на уровне \0-\3 %, а ошибка ложного обнаружения p'\ для изображений без искажений не превосходит 5 %.
Таким образом, для рассматриваемых систем ДЗЗ оптимальным является использование от \6 до 20 изображений (т.е. более 80 % от количества, необходимого для восстановления сцены), так как в этом случае обе ошибки p\ и p2 сохраняются на уровне \0 %.
Как видно из рис. 5, критическим количеством изображений является 4, так как в этом случае наблюдается резкий рост всех ошибок, за исключением p2 для Г < \000. Данный факт объясняется недостаточной длиной последовательности данных для анализа алгоритмом LOF при малом количестве изображений и малом размере суперпикселей. Увеличение статистики возможно за счёт увеличения размера суперпикселей, например при Г < \000. Однако это позволяет скомпенсировать только ошибку пропуска искажений, что видно на рис. 56, но не позволяет снизить ошибку ложного обнаружения.
Зависимость качества работы алгоритма от количества изображений с искажениями в составе сцены и без в наборе данных Другим важным фактором, влияющим на качество оценки искажений в составе сцены, является отношение 6 количества изображений с изменениями Н\ к общему количеству изображений Н:
6 = Е\/Н. (27)
Зависимость качества работы алгоритма от 6 оценивалась по выборкам изображений, полученным из базового набора путём исключения из него изображений без искажений. На рис. 6 приведены результаты эксперимента для 6 от 0,5 до 1, причём 6 = 1 соответствовало случаю, когда все изображения в наборе были с искажениями в составе сцены.
0,28 0,26 0,24 0,22 0,20 0,18 0,16 0,14 0,12 0,10 0,08
Р1 -о 500
о\ 1000 1500 2000
А \\
Л * А
\ «
г <ч
"Ч
Ь*—й
0,7 0,6, 0,5 0,4 0,3 0,2
Р2 -о— 500 -х— 1000 ■■■*■.......1500 -а-- 2000
] \
> \
1 \ \ \ \'\
\\ V» и
Г™1 •тагагг!
а)
8 10
12 14 16 18 20 Количество снимков
б)
8 10
12 14 16 18 20 Количество снимков
Р'1
\
-о— 500 — 1000
■■■*■.......1500
-о— 2000
в)
12 14 16 18 20 Количество снимков
Рис. 5. Зависимость качества работы алгоритма от общего количества изображений в наборе данных при различных значениях параметра Г: средняя ошибка ложного обнаружения на изображениях с искажениями в составе сцены (а), средняя ошибка пропуска искажений (б), средняя ошибка ложного обнаружения на изображениях без искажений в составе сцены (в)
0,08 0,07 0,06 0,05 0,04
ч -о— 1000 -х— 1500 ■■■*■.......2000
*
•ь
Vй
0,20
О,
а)
5 0,6 0,7 0,8 0,9 1,0 Доля изображений с искажениями
в составе сцены б)
0,15
0,10
0,6 0,7 0,8 0,9 1,0 Доля изображений с искажениями в составе сцены
Р'1 0,020
0,015 0,010 0,005 О
-О— 1000 -х— 1500 ■■■*.......2000 // Г
/ / / Р
Г% V \ /о' / у ✓
Ъ-СУ гУ Ч ч \ / у
в)
0,5 0,6 0,7 0,8 0,9
Доля изображений с искажениями в составе сцены
Рис. 6. Зависимость качества работы алгоритма от величины в: средняя ошибка ложного обнаружения на изображениях с искажениями в составе сцены (а), средняя ошибка пропуска искажений (б), средняя ошибка ложного обнаружения
на изображениях без искажений в составе сцены (в)
Из рис. 6 можно заключить, что Г практически не влияет на качество работы алгоритма при различном соотношении количества изображений с искажениями и без искажений. При любых значениях Г ошибка ложного обнаружения р1 снижается в среднем на 2% с ростом 6, а ошибка пропуска искажений р2 растёт в 2,5 раза. Ошибка ложного обнаружения на изображениях без изменений р'1 составляет не более 3 %.
Рост ошибки пропуска искажений связан с наличием в анализируемом наборе изображений, маски искажений для которых имеют общие части. То есть при малом количестве изображений без искажений пересекающиеся на различных снимках участки масок искажений могут быть не обнаружены.
Можно заключить, что для качественного обнаружения искажений в составе сцены предпочтительно, чтобы на каждое изображение с искажениями имелось одно изображение без искажений (6 = 0,5) или использовать изображения, участки с искажениями для которых не пересекаются. Очевидно, что если изображений без искажений будет больше (6 < 0,5), то ошибка р2 может быть значительно сокращена.
Другим вариантом снижения ошибки р2 при 6 > 0,5 является увеличение параметра у, отвечающего за порог аномальности для суперпикселей, частично занятых областью искажений. Для анализа влияния у на ошибки р1 и р2 был проведён эксперимент с набором изображений при 6 = 1 и рекомендуемых параметрах алгоритма. Результаты эксперимента представлены в табл. 4.
Табл. 4. Значения ошибок ложного обнаружения и пропуска искажений в составе сцены в зависимости от параметра у при 6 = 1 и рекомендуемых параметрах алгоритма
у р1 р2
0,1 0,071 0,233
0,15 0,087 0,195
0,05 0,065 0,235
0,2 0,112 0,160
0,25 0,143 0,126
Видно, что увеличение у с 0,1 до 0,25 позволяет в два раза снизить ошибку р2, однако при этом ошибка р1 также увеличивается в 2 раза. Тем не менее при остальных рекомендуемых параметрах алгоритма значение р1 не превосходит у. В практических случаях у следует выбирать как компромиссное значение между ошибками р1 и р2.
Зависимость качества обнаружения искажений в составе сцены от параметров пространственной и спектральной дискретизации изображений
Поскольку базовый набор изображений содержал изображения для двух систем ДЗЗ с различными параметрами дискретизации, то особый интерес представляло исследование влияния совместной обработки данных для двух систем ДЗЗ по отношению к результатам независимой обработки изображений для каждой из систем.
Для выявления влияния отличий в параметрах дискретизации изображений на результаты работы
алгоритма при рекомендуемых параметрах исходный набор из 20 изображений был разделён на два под-набора, содержащих 4 изображения для системы 1 и 16 изображений для системы 2. Соотношение изображений с искажениями и без искажений в обоих поднаборах было равно 1:1.
Результаты эксперимента приведены в табл. 5. Для удобства сравнения в табл. 5 в случае общего набора данных приведены два значения ошибок по изображениям системы 1 и по изображениям системы 2. Значения ошибок по всему набору данных для случая совместной обработки данных различных систем составили р\ = 0,088, р2 = 0,102, р'1 = 0,014.
Табл. 5. Оценка качества работы алгоритма на наборах данных с различными параметрами спектральной и пространственной дискретизации при рекомендуемых параметрах алгоритма
Обработанный набор данных Набор данных, использованный для расчёта ошибок р1 р2 р'1
Общий для двух систем 4 изображения системы 1 0,100 0,140 0,053
16 изображений системы 2 0,084 0,092 0,005
Система 1 4 изображения 0,098 0,455 0,054
Система 2 16 изображений 0,100 0,100 0,013
Из табл. 5 видно, что в случае системы 2 ошибки р1 и р2 при совместной обработке с изображениями системы 1 и при независимой обработке только изображений системы 2 близки. Причина близких значений ошибок - достаточное количество изображений для системы 2 (80 % от общего количества изображений в наборе). Тем не менее даже совместная обработка позволяет снизить для системы 2 значения всех ошибок, включая р'1, в среднем на 1-2 %. Для системы 1 с меньшим количеством изображений (4 изображения) совместная обработка позволила значительно снизить ошибку пропуска искажений, при этом ошибки р1 и р'1 остались практически на том же уровне.
Результаты эксперимента говорят о том, что совместная обработка данных с различными параметрами дискретизации позволяет сократить в среднем ошибки оценки искажений в составе сцены.
Оптимизация параметров алгоритма при малом количестве изображений
Дополнительно для случая малого количества изображений, а именно для изображений системы 1, был проведён эксперимент по подбору параметров алгоритма, обеспечивающих приемлемые (порядка 10%) ошибки пропуска искажений в составе сцены и ошибки ложного обнаружения.
Как было отмечено выше, уменьшение ошибки р2 требует уменьшения количества суперпикселей Г и увеличения у. Однако обе меры приведут к увеличению ошибки р1. Анализ участков изображений с ошибками ложного обнаружения р1 показал, что на малом количестве изображений (4 изображения) данные ошибки возникают для суперпикселей, частично
пересекающих области искажений в составе сцены. Для контроля ошибки ложного обнаружения в этом случае в алгоритм был введён параметр ю, задающий нижний порог доли пропускаемых аномальных кластеров в суперпикселе. Тогда возникающие лишние аномальные значения могут быть исключены из результата.
В табл. 6 для набора изображений системы 1 (4 изображения) и при различных значениях Г, у и ю приведены результаты оценки ошибок р1, р2 и р'1. Лучшие комбинации отмечены жирным шрифтом. Результаты позволяют заключить, что на малом количестве изображений путём подбора параметров Г, у и ю можно обеспечить значения всех ошибок на уровне 10-12 %.
Табл. 6. Оценка качества работы алгоритма при различных параметрах для малого количества изображений (4 изображения системы 1)
Г У ю p1 p2 p'1
2000 0,1 0,00 0,098 0,455 0,054
1500 0,1 0,00 0,113 0,367 0,085
1000 0,1 0,00 0,128 0,110 0,110
1000 0,1 0,25 0,113 0,122 0,098
500 0,1 0,25 0,111 0,147 0,068
500 0,2 0,25 0,113 0,098 0,071
500 0,2 0,50 0,109 0,128 0,067
500 0,3 0,50 0,107 0,136 0,067
Зависимость качества обнаружения искажений в составе сцены от доли пикселей с искажениями в составе сцены на каждом изображении
Очевидно, что при существенной степени пересечения областей с искажениями в составе сцены на анализируемых изображениях алгоритм не сможет распознать как искажения участки, подверженные одному и тому же роду искажений на подавляющем количестве изображений.
Для определения допустимой доли искажений л на каждом изображении при базовой конфигурации набора данных и рекомендуемых параметрах алгоритма был проведён эксперимент, результаты которого отражены в табл. 7.
Табл. 7. Оценка качества работы алгоритма при различной доле пикселей с искажениями в составе сцены
л p1 p2 p'1
10 % 0,087 0,100 0,014
20 % 0,069 0,190 0,002
30 % 0,058 0,288 0,000
40 % 0,059 0,335 0,000
50 % 0,103 0,408 0,000
Из табл. 7 видно, что при л = 20 ошибка пропуска искажений равна 19 % при рекомендуемом наборе параметров. Так как наибольшее влияние на ошибку р2 оказывает параметр у, был проведён дополнительный эксперимент для отыскания оптимальных значений у при прочих рекомендуемых значениях параметров с целью снижения ошибки р2 до 10 % при сохранении ошибки р1 на уровне 10 % для различных значений л на изображениях с искажениями. Результаты данного эксперимента представлены в табл. 8.
Табл. 8. Оптимальные значения у и показатели качества работы алгоритма при различных значениях л
У л p1 p2 p'1
0,1 10% 0,087 0,100 0,014
0,15 20% 0,106 0,117 0,009
0,2 30% 0,136 0,131 0,011
0,25 40% 0,169 0,137 0,009
0,3 50% 0,151 0,26 0,005
Результаты, приведённые в табл. 8, позволяют заключить, что при соответствующем подборе параметра у можно добиться компромиссного решения по ошибкам ложного обнаружения и пропуска искажений. Однако обеспечить их соответствие уровню 10-\3 % удаётся только для доли искажений в составе сцены не выше, чем 30 % на всех изображениях. Дальнейшее сокращение обеих ошибок возможно только за счёт увеличения количества изображений без искажений в анализируемом наборе данных.
Резюме экспериментальных исследований
Проведённые экспериментальные исследования позволяют заключить, что разработанный алгоритм способен оценивать маски искажений в составе сцены с ошибкой ложного обнаружения и пропуска искажений порядка 10-13 % в среднем. При этом доля искажений в составе сцены на каждом изображении не должна превышать 30 %. Алгоритм демонстрирует указанную точность обнаружения искажений для наборов данных с долей изображений с искажениями в составе сцены порядка 50 %. Увеличение доли неискажённых изображений в наборе приводит к увеличению точности оценки масок искажений.
Эксперименты по анализу совместного использования данных различного разрешения для оценки масок искажений показали преимущество объединённого набора данных с большим количеством изображений перед поднаборами с одинаковыми параметрами спектрального и пространственного разрешения, что свидетельствует о достижении положительного эффекта при комплексировании данных различного разрешения.
Заключение
В статье предложен алгоритм обнаружения случайных искажений в составе сцены по серии разновременных изображений ДЗЗ. Детектируемые искажения представляют собой локальные существенные изменения яркости наблюдаемой сцены, не связанные с изменением состояния наблюдаемых объектов, но проистекающие из изменения условий наблюдений территории в момент съёмки. Основу алгоритма составляют методы суперпиксельной сегментации и обнаружения аномалий. Первые позволяют учесть пространственную корреляцию пикселей изображений в серии, вторые - выделить спектрально-временные неоднородности в полученных локальных областях. В рекомендуемой реализации для получения разбиения на суперпиксели используется метод SLIC, а для обнаружения аномалий - метод LOF.
Экспериментальные исследования алгоритма на серии модельных изображений ДЗЗ с искажениями в
составе сцены показали, что при достаточном количестве данных для анализа алгоритм позволяет получить маски искажений в составе сцены с ошибкой ложного обнаружения и пропуска искажений порядка 10-13 % в среднем. При этом ошибка ложного обнаружения по изображениям без искажений не превышает 5 %. Данные показатели могут быть обеспечены для набора данных с общим числом изображений более 16 при доле пикселей с искажениями в составе сцены не более 30 % и отношении количества снимков без искажений к количеству снимков с искажениями 1:1.
Достоинствами алгоритма являются: универсальность относительно характеристик спектрального и пространственного разрешения анализируемых данных, отсутствие необходимости в привлечении априорных данных в виде метаданных снимков или данных метеорологических наблюдений, а также в оценке масок искажений в составе сцены с пространственным разрешением более высоким, нежели исходные изображения.
Предложенный алгоритм планируется применить в задаче восстановления изображения с повышением разрешения по серии разновременных изображений в качестве предварительной оценки релевантных для восстановления участков изображений.
Благодарности Работа выполнена при поддержке фонда РФФИ проекты № 18-07-00748 А, 16-29-09494 офи_м, а также содержит результаты проекта «Создание геоинформационного хаба Больших Данных», выполняемого в рамках реализации Программы Центра компетенций Национальной технологической инициативы «Центр хранения и анализа больших данных», поддерживаемого Министерством науки и высшего образования Российской Федерации по Договору МГУ имени М.В. Ломоносова с Фондом поддержки проектов Национальной технологической инициативы от 11.12.2018 № 13/1251/2018.
Литература
1. Pasetto, D. Integration of satellite remote sensing data in ecosystem modelling at local scales: Practices and trends / D. Pasetto, S. Arenas-Castro, J. Bustamante, R. Casagrandi, N. Chrysoulakis, A.F. Cord, A. Dittrich, C. Domingo-Marimon, G. El Serafy, A. Karnieli, G.A. Kordelas, I. Manakos, L. Mari, A. Monteiro, E. Palazzi, D. Poursanidis, A. Rinaldo, S. Rinaldo, S. Terzago, A. Ziemba, G. Ziv, G.A. Kordelas // Methods in Ecology and Evolution. - 2018. - Vol. 9, Issue 8. - P. 1810-1821.
2. Аншаков, Г.П. Комплексирование гиперспектральных и мультиспектральных данных КА «Ресурс-П» для повышения их информативности / Г.П. Аншаков, А.В. Ращупкин, Ю.Н. Журавель // Компьютерная оптика. - 2015. - Т. 39, № 1. - С. 77-82. - DOI: 10.18287/0134-2452-2015-39-1-77-82.
3. Белов, А.М. Спектральное и пространственное сверхразрешение при комплексировании данных ДЗЗ различных источников / А.М. Белов, А.Ю. Денисова // Компьютерная оптика. - 2018. - Т. 42, № 5. - С. 855-863. -DOI: 10.18287/2412-6179-2018-42-5-855-863.
4. Денисова, А.Ю. Алгоритмы анализа линейной спектральной смеси на гиперспектральных изображениях с использованием картографической основы / А.Ю. Денисова, В.В. Мясников // Компьютерная оптика. - 2014. - Т. 38, № 2. - С. 297-303.
5. Sun, L. A cloud detection algorithm-generating method for remote sensing data at visible to short-wave infrared wavelengths / L. Sun, X. Mi, J. Wei, J. Wang, X. Tian, H. Yu, P. Gan // ISPRS Journal of Photogrammetry and Remote Sensing. - 2017. - Vol. 124. - P. 70-88.
6. Thompson, D.R. Rapid spectral cloud screening onboard aircraft and spacecraft / D.R. Thompson, R.O. Green, D. Keymeulen, S.K. Lundeen, Y. Mouradi, D.C. Nunes, R. Castaño, S.A. Chien // IEEE Transactions on Geoscience and Remote Sensing. - 2014. - Vol. 52, Issue 11. - P. 67796792.
7. Zhu, X. An automatic method for screening clouds and cloud shadows in optical satellite image time series in cloudy regions / X. Zhu, E.H. Helmer // Remote Sensing of Environment. - 2018. - Vol. 214. - P. 135-153.
8. Li, J. High-spatial-resolution surface and cloud-type classification from MODIS multispectral band measurements / J. Li, W.P. Menzel, Z. Yang, R.A. Frey, S.A. Ackerman // Journal of Applied Meteorology. - 2003. - Vol. 42, Issue 2.
- P. 204-226.
9. Luo, S. Shadow removal based on clustering correction of illumination field for urban aerial remote sensing images / S. Luo, H. Li, H. Shen // 2017 IEEE International Conference on Image Processing (ICIP). - 2017. - P. 485-489.
10. Mostafa, Y. A review on various shadow detection and compensation techniques in remote sensing images / Y. Mostafa // Canadian Journal of Remote Sensing. - 2017.
- Vol. 43, Issue 6. - P. 545-562.
11. Movia, A. Shadow detection and removal in RGB VHR images for land use unsupervised classification / A. Movia, A. Beinat, F. Crosilla // ISPRS Journal of Photogrammetry and Remote Sensing. - 2016. - Vol. 119. - P. 485-495.
12. Champion, N. Automatic detection of clouds and shadows using high resolution satellite image time series / N. Champion // International Archives of the Photogramme-try, Remote Sensing & Spatial Information Sciences. -2016. - Vol. 41, Issue B3. - P. 475-479.
13. Breunig, M.M. LOF: identifying density-based local outliers / M.M. Breunig, H.P. Kriegel, R.T. Ng, J. Sander // ACM Sig-mod Record. - 2000. - Vol. 29, Issue 2. - P. 93-104.
14. Domingues, R. A comparative evaluation of outlier detection algorithms: Experiments and analyses / R. Domingues, M. Filippone, P. Michiardi, J. Zouaoui // Pattern Recognition. - 2018. - Vol. 74. - P. 406-421.
15. Zhao, H. Transformation from hyperspectral radiance data to data of other sensors based on spectral superresolution / H. Zhao, G. Jia, N. Li // IEEE Transactions on Geoscience and Remote Sensing. - 2010. - Vol. 48, Issue 11. - P. 39033912.
16. Farsiu, S. Fast and robust multiframe super resolution / S. Farsiu, M.D. Robinson, M. Elad, P. Milanfar // IEEE transactions on image processing. - 2004. - Vol. 13, Issue 10. - P. 1327-1344.
17. Student. The probable error of a mean / Student // Bio-metrika. - 1908. - Vol. 6, Issue 1. - P. 1-25.
18. Achanta, R. SLIC Superpixels / R. Achanta, A. Shaji, K. Smith, A. Lucchi, P. Fua, S. Süsstrunk // EPFL Technical Report. - 2010. - Vol. 149300. - P. 1-15.
19. Achanta, R. SLIC superpixels compared to state-of-the-art superpixel methods / R. Achanta, A. Shaji, K. Smith, A. Lucchi, P. Fua, S. Süsstrunk // IEEE Transactions on Pat-
tern Analysis and Machine Intelligence. - 2012. - Vol. 34, Issue 11. - P. 2274-2282.
20. Hartigan, J.A. Algorithm AS 136: A k-means clustering algorithm / J.A. Hartigan, M.A. Wong // Journal of the Royal Statistical Society. Series C (Applied Statistics). - 1979. - Vol. 28, Issue 1. - P. 100-108.
21. Margulis, D. Photoshop LAB color: The canyon conundrum and other adventures in the most powerful colorspace. / D. Margulis. - Peachpit Press, 2005.
22. Vane, G. The airborne visible/infrared imaging spectrometer (AVIRIS) / G. Vane, R.O. Green, T.G. Chrien, H.T. Enmark, E.G. Hansen, W.M. Porter // Remote Sensing of Environment. - 1993. - Vol. 44, Issues 2-3. - P. 127-143.
23. Bartos, M. Cloud and shadow detection in satellite imagery. Master thesis / M. Bartos // Prague: Czech Technical University in Prague Faculty of Electrical Engineering, 2017. - P. 1-57.
24. Xu, M. Cloud removal based on sparse representation via multitemporal dictionary learning / M. Xu, X. Jia, M. Pickering, A.J. Plaza // IEEE Transactions on Geosci-
ence and Remote Sensing. - 2016. - Vol. 54, Issue 5. -P. 2998-3006.
25. Schowengerdt, R.A. Remote sensing: models and methods for image processing / R.A. Schowengerdt. - Elsevier, 2006. - 560 p.
26. Spot-7 Satellite Sensor [Электронный ресурс]. - URL: https: //www. satimagingcorp. com/ satellite-sensors/ spot-7/ (дата обращения 09.07.2019).
27. Российская группировка ДЗЗ по состоянию на (01.03.2015 г.) // Геоматика. - 2015. - T. 1. - C. 106-112.
28. Farsiu, S. Multiframe demosaicing and super-resolution of color image / S. Farsiu, M. Elad, P. Milanfar // IEEE Transactions on Image Processing. - 2006. - Vol. 15. - P. 141-159.
29. Matlab [Электронный ресурс]. - URL: https://www.mathworks.com/products/matlab.html (дата обращения 09.07.2019).
30. Janssens, J.H.M. Outlier detection with one-class classifiers from ML and KDD / J.H.M. Janssens, I. Flesch, E.O. Postma // Proceedings of the Eighth International Conference on Machine Learning and Applications. - 2009. - P. 147-153.
Сведения об авторах
Белов Александр Михайлович, 1980 года рождения. В 2003 году с отличием окончил Самарский государственный аэрокосмический университет имени академика С.П. Королёва (Самарский университет) по специальности «Прикладная математика и информатика». В 2007 году получил степень кандидата физико-математических наук. В настоящее время работает доцентом кафедры геоинформатики и компьютерной безопасности в Самарском университете. Область научных интересов: обработка изображений и геоинформационные системы. Автор 37 научных публикаций, из них 14 статей в научных журналах. Член Поволжского отделения Российской ассоциации распознавания образов и анализа изображений.
E-mail: bam.post@gmail.com .
Сведения об авторе Денисова Анна Юрьевна см. стр. 853 этого номера.
ГРНТИ: 28.21.15, 28.17.19, 89.57.35, 89.57.45 Поступила в редакцию 09 июля 2019 г. Окончательный вариант - 20 сентября 2019 г.
Scene distortion detection algorithm using multitemporal remote sensing images
A.M. Belov1, A. Y. Denisova1 1 Samara National Research University, Moskovskoye Shosse 34, 443086, Samara, Russia
Abstract
Multitemporal remote sensing images of a particular territory might include accidental scene distortions. Scene distortion is a significant local brightness change caused by the scene overlap with some opaque object or a natural phenomenon coincident with the moment of image capture, for example, clouds and shadows. The fact that different images of the scene are obtained at different instants of time makes the appearance, location and shape of scene distortions accidental. In this article we propose an algorithm for detecting accidental scene distortions using a dataset of multitemporal remote sensing images. The algorithm applies superpixel segmentation and anomaly detection methods to get binary images of scene distortion location for each image in the dataset. The algorithm is adapted to handle images with different spectral and spatial sampling parameters, which makes it more multipurpose than the existing solutions. The algorithm's quality was assessed using model images with scene distortions for two remote sensing systems. The experiments showed that the proposed algorithm with the optimal settings can reach a detection accuracy of about 90% and a false detection error of about 10%.
Keywords: accidental scene-distortions detection, remote sensing image fusion, super-pixel image segmentation, anomaly detection.
Citation: Belov AM, Denisova AY. Scene distortion detection algorithm using multitemporal remote sensing images. Computer Optics 2019; 43(5): 869-885. DOI: 10.18287/2412-6179-201943-5-869-885.
Acknowledgements: The work was partly funded by the Russian Foundation for Basic Research under RFBR grants ## 18-07-00748 a, 16-29-09494 ofi_m and under the project "Creation of a Geographic Information Hub of Big Data", carried out as part of the Competence Center Program of the National Technological Initiative "Big Data Storage and Analysis Center", supported by the Ministry of Science and Higher Education of the Russian Federation under an agreement between M.V. Lomonosov Moscow State University and the Project Support Foundation of the National Technology Initiative, dated December 11, 2018 No. 13/1251/2018.
References
[1] Pasetto D, Arenas - Castro S, Bustamante J, Casagrandi R, Chrysoulakis N, Cord AF, Dittrich A, Domingo-Marimon C, El Serafy G, Karnieli A, Kordelas GA, Manakos I, Mari L, Monteiro A, Palazzi E, Poursanidis D, Rinaldo A, Rinaldo S, Terzago S, Ziemba A, Ziv G, Kordelas GA. Integration of satellite remote sensing data in ecosystem modelling at local scales: Practices and trends. Methods in Ecology and Evolution 2018; 9 : 1810-21.
[2] Anshakov GP, Raschupkin AV, Zhuravel YN. Hyperspectral and multispectral Resurs-P data fusion for increase of their informational content. Computer optics 2015; 39(1): 77-82.
[3] Belov AM, Denisova AY. Spectral and spatial superresolution method for Earth remote sensing image fusion. Computer Optics 2018; 42(5): 855-863. - DOI: 10.18287/2412-6179- 2018-42-5-855-863.
[4] Denisova AYu, Myasnikov VV. Algorithms of linear spectral mixture analysis for hyperspectral images using base map. Computer Optics 2014; 38(2): 297-303.
[5] Sun L, Mi X, Wei J, Wang J, Tian X, Yu H, Gan P. A cloud detection algorithm-generating method for remote sensing data at visible to short-wave infrared wavelengths. ISPRS journal of photogrammetry and remote sensing 2017; 124: 70-88.
[6] Thompson DR, Green RO, Keymeulen D, Lundeen SK, Mouradi Y, Nunes DC, Castaño R, Chien SA. Rapid spectral cloud screening onboard aircraft and spacecraft. IEEE Transactions on Geoscience and Remote Sensing 2014; 52: 6779-92.
[7] Zhu X, Helmer EH. An automatic method for screening clouds and cloud shadows in optical satellite image time series in cloudy regions. Remote Sensing of Environment 2018; 214: 135-53.
[8] Li J, Menzel WP, Yang Z, Frey RA, Ackerman SA. High-spatial-resolution surface and cloud-type classification from MODIS multispectral band measurements. Journal of Applied Meteorology 2003; 42: 204-26.
[9] Luo S, Li H, Shen H. Shadow removal based on clustering correction of illumination field for urban aerial remote sensing images. 2017 IEEE International Conference on Image Processing (ICIP) 2017; 485-9.
[10] Mostafa Y. A review on various shadow detection and compensation techniques in remote sensing images. Canadian Journal of Remote Sensing 2017; 43: 545-62.
[11] Movia A, Beinat A, Crosilla F. Shadow detection and removal in RGB VHR images for land use unsupervised classification. ISPRS Journal of Photogrammetry and Remote Sensing 2016; 119: 485-95.
[12] Champion N. Automatic detection of clouds and shadows using high resolution satellite image time series. International Archives of the Photogrammetry, Remote Sensing & Spatial Information Sciences 2016; 41(B3): 475-9.
[13] Breunig MM, Kriegel HP, Ng RT, Sander J. LOF: identifying density-based local outliers. ACM sigmod record 2000; 29(2): 93-104.
[14] Domingues R, Filippone M, Michiardi P, Zouaoui J. A comparative evaluation of outlier detection algorithms: Experiments and analyses. Pattern Recognition 2018; 74: 406-21.
[15] Zhao H, Jia G, Li N. Transformation from hyperspectral radiance data to data of other sensors based on spectral superresolution. IEEE Trans. Geosci. Remote Sens. 2010; 48(11): 3903-12.
[16] Farsiu S, Robinson MD, Elad M, Milanfar P. Fast and robust multiframe super resolution. IEEE transactions on image processing 2004; 13(10): 1327-44.
[17] Student. The probable error of a mean. Biometrika 1908; 1-25.
[18] Achanta R, Shaji A, Smith K, Lucchi A, Fua P, Susstrunk S. SLIC Superpixels. EPFL Technical Report 2010; 149300: 1-15.
[19] Achanta R, Shaji A, Smith K, Lucchi A, Fua P, Susstrunk S. SLIC superpixels compared to state-of-the-art superpixel methods. IEEE transactions on pattern analysis and machine intelligence 2012; 34(11): 2274-82.
[20] Hartigan JA, Wong MA. Algorithm AS 136: A k-means clustering algorithm. Journal of the Royal Statistical Society. Series C (Applied Statistics) 1979; 28(1): 100-8.
[21] Margulis D. Photoshop LAB color: The canyon conundrum and other adventures in the most powerful color-space. Peachpit Press; 2005.
[22] Vane G, Green RO, Chrien TG, Enmark HT, Hansen EG, Porter WM. The airborne visible/infrared imaging spectrometer (AVIRIS). Remote sensing of environment 1993; 44(2-3) : 127-43.
[23] Bartos M. Cloud and Shadow Detection in Satellite Imagery. Master Thesis. Czech Technical University in Prague Faculty of Electrical Engineering. Prague, 2017; 1-57.
[24] Xu M, Jia X, Pickering M, Plaza AJ. Cloud removal based on sparse representation via multitemporal dictionary learning. IEEE Transactions on Geoscience and Remote Sensing 2016; 54(5): 2998-3006.
[25] Schowengerdt RA. Remote sensing: models and methods for image processing. Elsevier; 2006.
[26] Spot-7 Satellite Sensor. Source: (https://wwwsatimaging-corp.com/satellite-sensors/spot-7/).
[27] Russian remote sensing constellation by date (01.03.2015). Geomatics 2015; 1: 106-112.
[28] Farsiu S, Elad M, Milanfar P. Multiframe demosaicing and super-resolution of color image. IEEE Transactions on Image Processing 2006; 15: 141-159.
[29] Matlab. Source: (https://www.mathworks.com/pro-ducts/matlab.html) .
[30] Janssens JHM, Flesch I, Postma EO. Outlier detection with one-class classifiers from ML and KDD. Proceedings of the Eighth International Conference on Machine Learning and Applications 2009; 147-53.
Author's information
Aleksandr Mikhailovich Belov (b. 1980) graduated from S.P. Korolyov Samara State Aerospace University (Samara University), majoring in Applied Mathematics and Informatics in 2003. He received candidate's degree in Physics and Mathematics in 2007. Currently he is holding position of associate professor of Geoinformatics and Computer Se-qurity department at Samara University. His research interests are currently focused on image processing and geoinfor-mation systems. He is author of 37 publications, including 14 papers. Member of the Russian Pattern Recognition and Image Processing Association. E-mail: bam.post@gmail. com.
The information about author Anna Yurievna Denisova you can find on page 585 of this isue.
Received July 09, 2019. The final version - September 20, 2019.
