ИССЛЕДОВАНИЕ ЗЕМЛИ ИЗ КОСМОСА, 2008, № 5, с. 37-41
МЕТОДЫ И СРЕДСТВА ОБРАБОТКИ ^^^^^^
И ИНТЕРПРЕТАЦИИ КОСМИЧЕСКОЙ ИНФОРМАЦИИ
УДК 528.88(15)
ОБНАРУЖЕНИЕ ИЗМЕНЕНИЙ НА ЗЕМНОЙ ПОВЕРХНОСТИ ПО РАЗНОВРЕМЕННЫМ МУЛЬТИСПЕКТРАЛЬНЫМ ИЗОБРАЖЕНИЯМ
© 2008 г. В. Н. Сагалович*, Э. Я. Фальков, Т. И. Царева
ФГУП "Государственный научно-исследовательский институт авиационных систем", Москва *Тел.: (499) 157-94-33; e-mail: sagalovich@gosniias.ru Поступила в редакцию 05.02.2008 г.
Предложен метод обнаружения изменений на местности, основанный на пороговой обработке разности нормированных разновременных изображений. Он инвариантен относительно линейных преобразований яркости анализируемых изображений и поэтому нечувствителен к изменениям калибровочных зависимостей сенсоров и к схемам радиометрической и атмосферной коррекции, описываемым линейными зависимостями. Для оценки эффективности метода использовались как результаты нахождения числа измененных пикселов в сопоставляемых синтезированных тестовых изображениях с гауссовыми шумами различной интенсивности, так и расчеты по разновременным цветным снимкам сканера TM ИСЗ Landsat площади повреждения, вызванной пожаром.
ВВЕДЕНИЕ
Существующие в настоящее время приборы и методы анализа данных дистанционного зондирования позволяют регулярно получать оценки характеристик земной поверхности и обнаруживать их изменения. Своевременное и точное определение этих изменений необходимо для понимания природы наблюдаемых явлений и для принятия обоснованных решений. Предложено множество подходов к исследованию этой проблемы, выбор которых зависит, в частности, от специфики изучаемой проблемы и не является однозначным [1].
Важным частным случаем является обнаружение изменений при обработке мультиспектраль-ных изображений одной и той же географической площади, полученных в моменты времени и ?2. Учитывая большой объем данных, собираемый сенсорами на борту летающих платформ, целесообразно предусмотреть этап автоматизированного обнаружения самого факта изменения на местности с последующим выявлением природы этого изменения компетентными специалистами. При этом необходимо учитывать, что к моменту ?2 могут быть другими и, вообще говоря, неизвестными характеристики атмосферы, угол Солнца, калибровка сенсора.
В настоящей статье разработан метод обнаружения изменений, инвариантный относительно линейных преобразований яркости анализируемых изображений и поэтому нечувствительный к изменениям калибровочных зависимостей сенсоров и к схемам радиометрической и атмосферной коррекции, описываемым линейными зависимостями. Инвариантность обеспечивается специальной нормировкой исходных изображений.
Метод основан на пороговой обработке разностей нормированных разновременных изображений в каждом канале наблюдения с последующим подсчетом суммарного числа измененных пикселов. Эффективность метода оценивается на примерах обнаружения изменений в синтезированных тестовых изображениях со случайными шумами различной интенсивности.
Следует отметить, что важная для практических применений инвариантность к линейным преобразованиям обеспечена в известном методе анализа разновременных мультиспектральных изображений, основанном на использовании канонического корреляционного анализа (multivariate alteration detection — MAD) [2]. Однако в этом методе отсутствует однозначное решающее правило, позволяющее провести сегментацию изображения и выделить измененные пикселы.
В статье сопоставляются результаты, полученные предложенным методом и методом MAD при определении по разновременным цветным изображениям сенсора TM ИСЗ Landsat площади, поврежденной пожаром.
РАСЧЕТНЫЕ СООТНОШЕНИЯ
Пусть Х1 и Х2 — монохроматические изображения размером I х К пикселов одной географической площади, полученные в моменты времени и ?2. На изображении Х2 ищутся те пикселы, которые можно считать изменившимися по сравнению с изображением Х1. Обозначим
= (X - теап(Х1))/ст1, *2 = (X? - теап(Х>))/ст2,
38
САГАЛОВИЧ и др.
где шеап(Х1), шеап(Х2) — средние значения яркости для изображений; сть а2 — соответствующие среднеквадратические отклонения.
Очевидно, что значения х1 и х2 инвариантны относительно линейных преобразований Х1 и Х2. При этом их средние значения шеап(х1) = = шеап(х2) = 0, дисперсия уаг(х1) = уаг(х2) = 1.
Введем переменную у, пробегающую значения от 0 до 255
y = (|x2 — xi| — min) х 255/(max — min),
(1)
g = о
;= T + i
Gi(T) = |Pi(T)1 X hg[g- И1 (T)]21 ,
255
J( T) = 1 + 2 [Pi(T) lnGl + P2(T) ln02(T)] -- 2[ Pi( T) ln Pi( T) + P2 (T) ln P2 (T) ].
(2)
Метод Капура Обозначим H(T) = — к. 1пй„. Искомое T
= 0 6 а
вычисляется из условия максимума функционала ¥(7):
¥( T) = ln [ Pi ( T) P2 ( T) ] + H( T) /Pi ( T) + + [H(255) - H( T)]/P2( T).
(3)
где min, max — минимальное и максимальное значения модуля разности x2 — x1.
При пороговой обработке [3] изображения y ищется подходящее значение порога T такое, что при y,k > T(i = 0, 1, ... I- 1; к = 0, 1, ... K- 1) пикселы соответствуют измененной части изображения. Предложены различные подходы к расчету T, однако накопленный опыт обработки изображений показал [3], что предпочтительными оказываются известные методы минимальной ошибки Киттлера [4] и максимальной энтропии Капура [5]. Поэтому в данной статье использовались именно эти методы.
Для выполнения расчетов необходимо построить гистограмму hg (g = 0, 1, ... 255) выборки значений яркости изображения у.
Метод Киттлера Введем обозначения:
T 255
Pi ( t = X ^ P2 ( T = X ^
g = 0 g = T + i
T 255
Hi(T) = Pi(T)-i Xhgg, T) = P2(T)-i X hgg,
После определения величины порога Т для конкретного канала мультиспектрального изображения можно найти количество измененных пикселов. Суммируя по всем каналам и учитывая каждый измененный пиксел только один раз, находим их общее количество, которое, в частности, может служить критерием при автоматизированном решении вопроса о передаче изображения на рассмотрение специалистам для дальнейшего анализа.
Анализ тестовых изображений
Для оценки эффективности различных алгоритмов сегментации предпочтительно использовать синтезированные тестовые изображения, поскольку они легко воспроизводятся в различных исследованиях и их характеристики просто контролировать [6]. В настоящей работе первое из разновременных изображений задавалось формулой
Imp = 112 + 32sin(p/«)sin(q/«),
G2( T) = Ы( T)-i X hg[g - ( T) ]2 j- .
^ g = T + i '
Искомое значение T вычисляется из условия минимума функционала J(T)
(4)
гдер = 0, 1, ... 255 и q = 0, 1, ... 255, а п характеризует размеры темных и светлых областей на изображении. Второе изображение отличалось от первого только тем, что в его центре размещался круг радиуса г с уровнем яркости 144, причем принималось, что в формуле (4) п = 2п/г. На рис. 1 представлены эти изображения (при г = 32.3).
При оценке точности расчетов по определению площади изменений 8с следует учитывать, что значение этой площади для тестовых изображений известно и составляет яг2. Ошибка в определении 8с оценивается по формуле Ег = Пг2 — ^с|/яг2 [6].
Были проведены расчеты Ег для разновременных тестовых изображений при наложении на них гауссовых шумов с нулевым средним и сред-неквадратическими отклонениями а = 0, 4, 8 и 16. Величины порогов Т находились методами Киттлера и Капура с помощью оптимизации функционалов (2) и (3). Кроме того, простым перебором находилось оптимальное значение Т, которое минимизировало ошибку Ег.
На рис. 2 представлены зависимости ошибок оценки площади изменений Ег от уровня шумов. Видно, что при наличии шума метод Капура менее эффективен, чем метод Киттлера, и оба эти метода не позволяют достигнуть результатов, полученных при оптимальном выборе значения по-
ОБНАРУЖЕНИЕ ИЗМЕНЕНИЙ НА ЗЕМНОЙ ПОВЕРХНОСТИ
39
&
ь Иг
и ■
^чци
б
Рис. 1. Синтезированные тестовые изображения: а — исходное изображение; б — изображение с измененной яркостью внутри центрального круга (г = 32.3).
рога Т В дальнейших расчетах порог Т будет находиться методом Киттлера.
Выявление изменений, вызванных пожаром, по разновременным цветным изображениям
С целью апробации предложенного метода применительно к натурным изображениям были обработаны цветные снимки, полученные с помощью сканера ТМ ИСЗ Ьапё8а1 для области вбли-
Ошибка
Ошибка
ст
Рис. 2. Ошибки определения площади изменений по разновременным тестовым изображениям с шумами различной интенсивности (a1 — оптимальное значение порога, a2 — метод Киттлера, a3 — метод Капура): а — задаваемая относительная площадь изменений 0.05; б — задаваемая относительная площадь изменений 0.15.
зи оз. Тахоа, шт. Невада, до пожара (8.05.1986 г.) и после пожара (8.05.1992 г), размещенные на Интернет-сайте центра данных EROS [7] (соответствующие полутоновые изображения в цветовой системе NTSC ([8], стр. 216) представлены на рис. 3). С помощью стандартных алгоритмов по этим снимкам были получены монохроматические изображения (размером 200 х 200 пикселов, площадью Sim = 4 х 104) в соответствующие моменты времени для красного, зеленого и синего каналов. Для каждого канала рассчитывалась разность нормированных изображений, которая преобразовывалась по формуле (1), и строилась гистограмма значений у. Далее из условия минимума функционала J(T) (2) находились значения порогов и площади Sc, занятые измененными пикселами. Было получено, что отношение SJSim равно для красного канала 0.14, для зеленого канала — 0.13 и для синего — 0.11.
При нахождении общей площади поражения пожаром Sf оказалось, что многие пикселы по-
40
САГАЛОВИЧ и др.
Рис. 3. Полутоновые изображения области вблизи оз. Тахоа, шт. Невада, соответствующие цветным снимкам, полученным с помощью сканера TM ИСЗ Landsat: а — снимок от 08.05.1986 г. (до пожара); б — снимок от 08.05.1992 г. (после пожара).
вторяются как измененные в двух или трех каналах. Эти пикселы учитывались только один раз. Расчеты показали, что S/Sim = 0.23.
Для нахождения ошибок Er использовались результаты обработки представленного в [7] дополнительного изображения, на котором выделена площадь поражения Sa, вызванного пожаром. Было
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.