Одне з головних завдань дистанційного зондування Землі – це моніторинг об’єктів та явищ на її поверхні. Регулярність космічної зйомки дозволяє вести постійне спостереження та оперативно виявляти й оцінювати зміни в природних та перетворених людиною ландшафтах. Це дуже корисна інформація, яка знаходить своє застосування в контролі землекористування, лісокористування, боротьбі з наслідками надзвичайних ситуацій та інших галузях. Регулярна космічна зйомка виконується вже декілька десятиріч. За цей час накопичено величезний архів космічних знімків. Це фактично літопис стану ландшафтів Землі. Він дозволяє виявити зміни через роки після того як вони сталися, тобто провести їх ретроспективний аналіз. Серед сучасних технологій обробки космічних знімків важливе місце займають технології автоматизованого виявлення та картографування змін за різночасовими космічними знімками. Такі технології є головним елементом системи дистанційного моніторингу. Цим постом ми плануємо відкрити цілу серію публікацій про технології дистанційного виявлення змін за допомогою космічних знімків.
Теорія
Космічні знімки – це не тільки графічне зображення земної поверхні, але також її кількісне описання. Кожен піксель знімку має кількісне значення – значення відбиття сонячної радіації від земної поверхні. Тому космічні знімки можна обробляти за допомогою математичних методів. Математичні операції – це одна з найпростіших технологій виявлення змін на космічних знімках. Ці операції є окремим випадком математичної обробки растрів взагалі. І ніяких специфічних рис не мають. Для виявлення змін використовують віднімання та ділення. Операцію віднімання растрів можна описати наступною формулою:
У цій формулі I1 та I2 – це растри двох дат (літера I від слова «image», тобто зображення). D – це зображення різниці, «difference» – різниця. Літери x та у нижньому індексі – це координати пікселя. Літера b y нижньому індексі використається для багатошарових растрів і позначає номер шару. C – це константа, яку додають щоб оминути негативних значень різниці. Але часто негативні значення не сприймаються як щось погане, тому зазвичай таку константу не додають.
З наведеної формули ми бачимо, що від значення пікселя за одну дату віднімається значення пікселя з такими самими координатами за іншу дату. В результаті зображення різниці має такі ж самі характеристики, як і основні растри. Якщо ці растри багатошарові та мають однакову кількість шарів, то зображення різниці теж буде багатошаровим і з такою ж кількістю шарів (рис. 1).
Рис. 1. Принцип віднімання растрів
Операцію ділення космічних знімків можна описати наступною формулою:
На відмінність від формули віднімання, для ділення результат позначений літерою Q (від слова «quotient»). Константа у знаменнику використовується щоб оминути ділення на нуль.
Растри з різною кількістю шарів теж можна віднімати та ділити. Але для космічних знімків це не має сенсу, бо знімки, що зроблені одним сенсором, мають однакову кількість каналів. А порівнювати знімки різних сенсорів – не найкраще рішення.
Коли ми виконуємо віднімання або ділення треба забезпечити коректність цих процедур. Це значить що виявлені зміни яскравості земної поверхні мусять відповідати реальним змінам земної поверхні. Для цього краще використати знімки, що зроблені одним сенсором. Тоді відмінності яскравості на знімках будуть пов’язані зі змінами земної поверхні, а не з різницею в технології зйомки. Якщо треба порівняти стан місцевості за різні роки, бажано щоб знімки були зроблені в той самий сезон. Тоді ми гарантовано виявимо багаторічні зміни, а не сезонну динаміку місцевості.
Після підбору знімків для порівняння необхідно ще виконати їх радіометричне калібрування та атмосферну корекцію. Завдяки цьому різниця в стані атмосфери не буде заважати виявляти реальні зміни на місцевості.
Практика
Практику пошуку змін місцевості за космічними знімками в програмі ENVI ми розберемо на прикладі виявлення вирубки лісів. На рисунку 2 показані два фрагменти знімків із супутника Landsat 5TM (номер сцени path/row 177/025). Вони охоплюють один і той же масив листяного лісу, що розташований у Вовчанському (північна частина) та Печенізькому (південна частина) районах Харківської області. Цей ліс протягнувся вздовж правого берега Печенізького водосховища, від Старого Салтова до Печеніг. Північна його частина має назву Хотімлянська дача, а південна – Печенізька дача.
Рис. 2. Фрагменти знімку Landsat 5 TM за 1996 (вгорі) та 2007 рік (внизу)
Перший знімок був знятий 7 травня 1996 року. Другий знімок був знятий через одинадцять років – 6 травня 2007 року. Тобто обидва знімки було отримано в один сезон, а фенофази розвитку рослинності мають бути майже однаковими. Тому там де в покриві лісу не було змін, значення яскравості пікселів будуть відрізнятися слабко. Звісно, змінився вік лісу. Але відповідна цьому зміна яскравості буде відносно невеликою та майже однаковою для всіх частин лісу. А там де трапилися значні зміни місцевості, зміни яскравості будуть великими. Такими змінами місцевості можуть бути лісові пожежі, вируб лісу, всихання лісу внаслідок поширення хвороб та шкідників лісу.
Візуально на знімках (рис. 2) ми бачимо, що здебільшого ліс в 1996 та 2007 роках виглядає однаково. Але для 2007 року впадають в очі чисельні вирубки, яких немає в 1996 роки. Наша мета – дешифрувати таких ділянок в автоматизованому режимі.
Виявляти зміни на території лісу будемо за допомогою операції віднімання знімків. Ділення в ENVI виконується подібним зразком, тому окремо його ми розглядати не будемо. Питання радіометричного калібрування та атмосферної корекції в рамках цього посту теж розглядатися не будуть. Це окремі та доволі великі теми.
Алгоритм, за яким ми будемо виявляти вирубки, показаний на рисунку 3. Він складається з чотирьох кроків:
1) Створення з каналів знімку усередненого (псевдопанхроматичного) зображення
2) Віднімання різночасових знімків
3) Створення карти класифікації за допомогою квантування
4) Конвертування карти класифікації в векторний формат
Перший крок вимагає пояснень. Він не є обов’язковим. Коли ми віднімаємо з одного знімку іншій, ми отримуємо різницеве зображення для кожного каналу. Далі можна окремо обробляти кожне різницеве зображення. А можна визначити канал, для якого різниця найкраще передає зміни місцевості. Але щоб спростити собі завдання робиться перший крок нашого алгоритму – з каналів багатоспектральних знімків створюються усередненні зображення. Такі зображення ще називають псевдопанхроматичними. Як і справжні панхроматичні знімки вони складаються лише з одного шару, візуалізуються в чорно-білій гамі та містять інформацію з широкого діапазону спектру. Але панхроматичні знімки дійсно зняти сенсором, чутливим до світла у широкому діапазоні спектра. А псевдопанхроматичні зображення зробленні через усереднення декількох зображень в кількох вузьких діапазонах спектру.
Рис. 3. Схема обробки даних
1. Створення псевдопанхроматичних знімків
1) Щоб створити усереднене зображення, в Тулбоксі оберіть команду Statistic → Sum Data Bands. З’явиться вікно Sum Data Parameters (рис. 4). Це інструмент, що вираховує статистичні характеристики набору каналів знімка для кожного пікселя та створює відповідні зображення – суму, середнє, стандартне відхилення каналів та інше.
2) В верхній половині вікна є список цих параметрів (Select Outputs Bands:). В нашому випадку в списку треба обрати лише Mean – середнє значення.
3) Останній крок – обрати спосіб збереження результату та натиснути ОК.
На всіх проміжних кроках обробки даних краще зберігати результати в тимчасову пам’ять (варіант Memory для параметра Output Result to). І лише кінцевий результат – карту класифікації – зберігати в файл. Потім всі непотрібні результати буде стерто з тимчасової пам’яті під час вимикання програми.
Рис. 4. Налаштування розрахунку усередненого зображення
Результати розрахунку усереднених (псевдопанхроматичних) зображень наведено на рисунку 5. Вирубки, що з’явилися до 2007 року, там видно так само добре, як і на багатошарових знімках (рис. 2).
Рис. 5. Усередненні зображення за 1996 та 2007 рік, що обрізані по межі лісу
2. Створення різницевого зображення
1) Щоб почати створення різницевого зображення, оберіть в Тулбоксі ENVI команду Band Ratio→Band Math. Після цього з’явиться вікно математичних операцій над растрами (рис. 6, ліворуч).
2) Головний елемент вікна Band Math – це рядок вводу Enter an expression:, який стоїть у середній частині вікна. В рядок вводять математичні вираження для обробки растрів. Щоб отримати растр різниці, треба ввести вирази: b2-b1. Літерою b позначають змінні (канали знімку), а цифрою – номер каналу. Порядок нумерування каналів значення не має.
Для ділення растрів треба ввести в рядок вводу Enter an expression вираз: b2/b1. Це єдина відмінність, за якою в ENVI відрізняються операції віднімання та ділення растрів.
3) Після того, як вираз введено, треба натиснути кнопку Add to List. Якщо в синтаксисі немає помилок, то вираз буде додано у перелік (Previous Band Math Expressions), що розташований у верхній частині вікна Band Math.
Під переліком виразів знаходяться чотири кнопки, призначені для керування цими виразами. Кнопка Save дозволяє зберегти вираз в файлі з розширенням exp. Відкрити вираз з такого файлу можна за допомогою кнопки Restore. Кнопка Clear дозволяє очистити весь перелік виразів разом, а кнопка Delete дозволяє видалити з нього окремий обраний вираз.
4) Якщо вираз було додано до переліку Previous Band Math Expressions, можна натиснути кнопку ОК. Після цього виникне вікно (Variables to Bands pairings), у якому треба встановити відповідність між змінними з виразу та растрами (рис. 6, праворуч).
5) У верхній частині вікна Variables to Bands pairings розташований перелік змінних, що використовуються у виразі (Variables used in expression). Нижче – перелік растрів (Available Bands list). Щоб зіставити їх, треба обрати лівою кнопкою миші в переліку спочатку змінну, а потім відповідний їй растр.
Рис. 6. Вікно завдання математичного виразу (ліворуч) та вікно зіставлення змінних з растрами (праворуч)
Зверніть увагу, на те, що на рисунку 6 немає налаштувань збереження вихідних даних. Це тому що ще не всі пари змінна-растр зіставлені. Коли це буде виконано, у вікні з’явиться опція вибору способу збереження (в тимчасову пам’ять або в файл).
В результаті отримуємо зображення віднімання та ділення знімків, вони наведені на рисунку 7. На них добре видні вируби. Вони виглядають як чорні ділянки. На їх території між 1996 та 2007 роками відбулося зменшення яскравості поверхні.
Рис. 7. Зображення різниці (ліворуч) та зображення частки (праворуч), що обрізані по межі лісу
3. Квантування різницевого зображення
1) Щоб почати квантування треба в Тулбоксі ENVI обрати команду Classification → Raster Color Slices. Після цього з’явиться вікно, у якому можемо обрати зображення для обробки. В нашому випадку це зображення різниці.
2) Після обрання зображення для обробки з’явиться вікно налаштувань процедури квантування – Edit Raster Color Slices (рис. 8,9). В ньому будуть вже встановлені класи за замовчуванням. Їх необхідно видалити, що робиться за допомогою кнопки 
. Після цього вікно на вигляд стане як на рисунку 8. Ліворуч буде пуста таблиця – перелік класів. Таблиця має три стовпчики: колір класу (Color), нижня межа класу (Slice Min), верхня межа класу (Slice Max). Праворуч у таблиці стоїть гістограма значень растру.
Рис. 8. Вікно налаштувань процедури квантування без завданих класів
3) Далі можна завдати свої класи. Вони додаються натисканням кнопки 
. Нам треба додати два класи. Перший клас – це будуть вирубки, другий клас – ліс, що не зазнав великих змін.
4) Коли класи додано, необхідно виправити межі класів, що позначені в таблиці. Вони задаються вручну, з клавіатури. Оптимальні значення треба визначити користувачу шляхом спроб. Після додавання класів та налаштування їх меж вікно на вигляд буде як на рисунку 9.
Рис. 9. Вікно налаштувань процедури квантування з класами, що задані користувачем
5) Далі можна налаштувати кольори класів. Для цього треба натиснути правою кнопкою миші по відповідній чарунці в стовпчику Color. З’явиться контекстне меню, у якому треба обрати команду Edit Color… Ця команда викликає вікно налаштування кольору (рис. 10). В ньому можна налаштувати колір за допомогою різних способів: – обрати зі списку, встановити за допомогою бігунців або вказати в кольоровій таблиці. У нашому прикладі для вирубів був встановлений червоний колір, а для інших частин лісу – синій.
Рис.10. Вікно налаштування кольору
6) Після того, як Ви натиснете кнопку ОК у вікні квантування, з’явиться результат квантування – карта класифікації. Вона буде збережена в тимчасовій пам’яті. Результат квантування для нашого прикладу наведений на рисунку 11. Виявлено вируби загальною площею 216 гектарів, що складає 2,3% від площі лісового масиву.
Рис. 11. Результат квантування (обрізаний по межі лісу)
7) Щоб зберегти карту класифікації у файл необхідно в таблиці змісту (Layer Manager) натиснути правою кнопкою миші на карту класифікації. З’виться контекстне меню, у якому можна вибрати команду Export Color Slices → Class Image… або Export Color Slices → Shapefile. Перша команда зберігає растр. Друга команда зберігає векторний шар у формати шейп-файлу.
Рис.12. Збереження результатів квантування
На цьому процес виявлення змін закінчується. Збереження в шейп-файл потрібно для подальшої обробки результатів в ГІС. Використання ГІС дозволяє створити якісні карти змін місцевості. Також можна провести поглиблений просторовий аналіз та описати закономірності розташування виявлених змін. А це вже крок на шляху до прогнозування розвитку території та пропонування заходів для запобігання небажаних змін.
