Автоматизований пошук змін на різночасових космічних знімках: метод діаграми розсіювання

Аналіз багатомірного простору спектральних ознак – це дуже корисний засіб для дешифрування багатоканальних космічних знімків. Його використання допомагає обрати оптимальний алгоритм класифікації з навчанням. Також можна безпосередньо в багатомірному просторі спектральних ознак виокремлювати класи об’єктів, а потім переносити їх на двомірний географічний простір космічного знімка. Ці обидва напрямки використання пов’язані з обробкою знімків однієї дати зйомки. Але простір спектральних ознак можна застосовувати також для обробки декількох знімків, що знято в різні дати. Це один зі способів автоматизованого пошуку змін за даними космічної зйомки. Ми розберемо його в цьому пості на найпростішому прикладі з двома різночасовими знімками.

 

Теорія

Щоб порівняти значення пікселів в одному каналі зйомки за дві різни дати, треба створити двомірний простір спектральних ознак, тобто двомірний графік, осі якого відповідають датам зйомки. Одна його вісь показує значення за першу дату, друга – за другу дату. Логічним буде, якщо вісь X (вісь абсцис) – це початкова дата, а вісь Y (вісь ординат) це кінцева дата. Трохи пояснемо цю логіку. Стан місцевості в кінці є результатом перетворень, що зазнав початковий її стан. Математично кажучи, фінальний стан (Y) є функцією від початкового стану (X): Y = f(X). А коли ми відображаємо графік функції, то відкладемо Y на вертикальній вісі, а X – на горизонтальній.

Всі можливі співвідношення початкових та фінальних значень в двомірному просторі виглядають як розсіяні точки. Тому такий графік називають також діаграмою розсіювання (так ми і будемо називати його далі). Інколи його називають скаттерограмою або скаттерплотом, від англійського scatter plot. Його схематичний приклад показаний на рисунку 1.

Рис. 1. Діаграма розсіювання

 

Значення, які між двома датами змінювалися слабо, створюють на діаграмі хмару значень, яка має вузьку та витягнуту форму (на рисунки 1 вона обведена зеленим пунктиром). Ця хмара витягнута вздовж лінії, яку можна описати формулою Y = kX+b. Коефіцієнти k та b показують відповідно кут нахилу лінії та її зміщення по вертикалі. В ідеалі,  коли немає змін місцевості k = 1 та b = 0. У такому разі лінія йде через початок координат і має кут нахилу 45º. Тобто для співпадаючих у просторі пікселів значення спектральної яскравості за обидві дати однакові.

Бувають випадки коли на місцевості змін немає, а зміни яскравості пікселів є. Це буває пов’язане з відмінностями в умовах зйомки, або з тим, що порівнюють дані від різних супутникових сенсорів. Якщо такі зміни яскравості пікселів однакові для всього знімка, то вони не заважають автоматизованому пошуку змін місцевості.  Пікселі, в яких не було змін, будуть створювати лінійно витягнуту хмару на діаграмі розсіювання. Але у такому разі вона уже буде мати кут нахилу, відмінний від 45º.

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

 

Практика

Практичні аспекти пошуку змін розберемо на прикладі, що ілюстрував один з минулих постів. Наше завдання – по знімках зі супутника Landsat виявити вирубки в лісовому масиві, що розташований на правому березі Печенізького водосховища.  Для цього ми проаналізуємо два знімки від 7 травня 1996 року та 6 травня 2007 року. Щоб діаграма розсіювання була більш наочною та для полегшення її інтерпретування, будемо використовувати знімки, які обрізані за межею лісу.

Рис. 2. Територія дослідження. Ліворуч – фрагмент знімку від 7.05.1996, праворуч – фрагмент знімку від 6.05.2007

 

Схема процесу обробки космічних знімків наведена на рисунку 3. Вона складається з п’яти етапів:

1) Створення із кожного знімку псевдопанхроматичного зображення за допомогою усереднення пікселів із різних каналів

2) Об’єднання двох різночасових зображень в один багатошаровий растр (пакетування растрових шарів)

3) Створення діаграми розсіювання та виокремлення на ній класу, що відповідає змінам місцевості

4) Експорт в область обробки класу, якій був виокремлений на діаграмі розсіювання

5) Експорт області обробки в векторний шар.

Перший та п’ятий етапи вже описувалися в минулому пості. Тому зараз ми не будемо їх докладно розбирати.

Рис. 3. Схема процесу обробки даних

 

Пакетування растрів

А) Програма ENVI дозволяє створювати двомірні діаграми розсіювання лише для растрових шарів, що належать одному файлу. Точніше це стосується сучасних її версій. У версіях старіше, ніж 5.0 можна було робити це також для растрових шарів з різних файлів. Треба було лише щоб вони мали однаковий розмір у пікселях та однакові дисплейні координати.

Під час обробки різночасових знімків ми спочатку маємо щонайменше два окремі знімки. Ми зробили з окремих знімків два окремих різницевих зображення. Далі їх необхідно об’єднати в одному багатошаровому растрі (пакеті растрових шарів). Для цього в Тулбоксі ENVI треба обрати команду Raster Manager → Layer Stacking. Далі отримаємо вікно параметрів пакетування (рис. 4).

Б) За допомогою кнопки Import File… треба вказати растри для пакетування. Після її натискання, з’явиться вікно, у якому можна обрати ці растри.

В) Далі треба налаштувати спосіб збереження результату (в тимчасову пам’ять або в постійний файл).

Г) Інші налаштування можна залишити такими, як є за замовчуванням.

Рис. 4. Вікно налаштування процесу пакетування растрових шарів в єдиний файл

 

Робота з діаграмою розсіювання

А) Щоб створити діаграму розсіювання треба в головному меню ENVI обрати команду Display → 2D Scatter Plot. Після цього з’явиться вікно двомірної діаграми розсіювання (2D Scatter Plot). Праворуч у цьому вікні знаходиться сама діаграма розсіювання, ліворуч – її налаштування та опції створення класів (рис. 5).

Рис. 5. Двомірна діаграма розсіювання

 

Б) Перш за все налаштуванню підлягають осі графіку. Розгортувані списки X Band та Y Band дозволяють завдати шари, дані з яких відображаються відповідно за горизонтальною та вертикальною віссю. В нашому прикладі (рис. 5) вісь X – це дані 1996 року, а вісь Y – це дані 2007 року.

В) Після налаштування осей треба вказати, дані з якої області знімку відображати на діаграмі. Є два варіанти. Можна використати дані зі всього знімку, а можна лише з тієї його частини, що зараз відображена у вікні ENVI. Якщо поставити галочку для параметра Viewable Area Only, ми оберемо другий варіант, якщо знімемо її – перший варіант.

Г) На діаграмі розсіювання окремою точкою відображається кожне співвідношення значень за першу та другу дату. Якісь співвідношення можуть траплятися в багатьох пікселях. Інші, навпаки, є тільки в нечисленних пікселях. Відобразити ці відмінності можна (щільність значень), якщо кількість пікселів для кожного співвідношення на діаграмі розсіювання позначити кольором. За замовчуванням у вікні 2D Scatter Plot щільність значень не відображається. Всі точки на графіку мають однаковий чорний колір. Якщо поставити галочку для параметра Density Slice, точки отримають колір. Точки з високою щільністю значень будуть мати зелений колір. Зі зменшенням щільності значень колір буде змінюватися на блакитний, потім на синій і для найменшої щільності він стане фіолетовим.

Рис. 6. Звичайна діаграма розсіювання (ліворуч) та діаграма розсіювання с позначеною кольорами щільністю значень

 

Д) Тепер розберемо створення класів (Class Functions). Опції створення класів окреслено чорною рамкою. Щоб виокремити клас у двомірному просторі, треба спочатку обрати у розгортуваному списку колір. Кожний клас, який ми виокремлюємо, мусить бути свого окремого кольору. За замовчуванням встановлений червоний колір.

Далі треба кліком лівої кнопки миші обводити частину хмари значень, яка відповідає пікселям зі змінами. Щоб закінчити цей процес, треба клацнути правою кнопкою миші.

Результат цього показаний на рисунку 7. Якщо клас обведено погано, то його можна прибрати натиснувши кнопку Clear, а потім обрати його знов. Після виокремлення класу змін можна його експортувати в область обробки. Це роблять за допомогою кнопки Export. Саме ця процедура дозволяє ідентифікувати розташування в географічному просторі точок, які ми виокремили.

Рис. 7. Двомірна діаграма розсіювання з позначеним червоним кольором класу вирубів

 

Робота з областю обробки та експорт її в векторний шар

Коли ми маємо область обробки, ми можемо накласти її поверх знімку. Приклад цього показаній на рисунку 8, де межі вирубів окреслені червоною лінією.  

Рис. 8. Виявлені вирубки у вигляді області обробки (позначена червоним), яка накладена згори космічного знімка. Ліворуч – знімок 1996 року, праворуч – знімок 2007 року

Для подальшого використання результатів дешифрування, наприклад, для створення карти змін за допомогою ГІС-програм, треба експортувати область обробки в векторний шар. Для цього треба в меню вікна областей обробки вибрати команду File → Export →Export to Shapefile… . Також результати дешифрування можна передати в ГІС в растровому форматі. Для цього область обробки треба конвертувати в карту класифікації. Ця процедура виконується командою Region of Interest → Classification Image from ROIs.