4. Короткий вступ до дистанційного зондування

4.1. Основні поняття

This chapter provides basic definitions about GIS and remote sensing.

4.1.1. Поняття ГІС

There are several definitions of GIS (Geographic Information Systems), which is not simply a program. In general, GIS are systems that allow for the use of geographic information (data have spatial coordinates). In particular, GIS allow for the view, query, calculation and analysis of spatial data, which are mainly distinguished in raster or vector data structures. Vector is made of objects that can be points, lines or polygons, and each object can have one or more attribute values; a raster is a grid (or image) where each cell has an attribute value (Fisher and Unwin, 2005). Several GIS applications use raster images that are derived from remote sensing.

4.1.2. Поняття дистанційного зондування

В загальному розумінні дистанційне зондування це «наука та технологія за допомогою яких характеристики об’єктів інтересу можуть ідентифікуватись, вимірюватись та аналізуватись за відсутності безпосереднього контакту» (JARS, 1993).

У вужчому розумінні дистанційне зондування це вимірювання енергії, що випромінюється від земної поверхні. Якщо джерелом енергії, що вимірються, є Сонце, то це пасивне дистанційне зондування і результатом такого вимиріювання може бути цифровий знімок (Richards and Jia, 2006). Якщо енергія, що вимірюється, випромінюється не Сонцем, а платформою сенсора, такою як радарні сесори, що працюють у мікрохвильовому діапазоні, то це активне дистанційне зондування (Richards and Jia, 2006).

Електромагнітний спектр це «система, що класифікує за довжиною хвилі всю енергію (від короткохвильової космічної до довгохвильової радіо), що гармонічно рухається з постійною швидкістю світла» (NASA, 2013). Пасивні сенсори вимірюють енергію з оптичних областей електромагнітного спектра: видимої, ближньої інфрачервоної (ІЧ), короткохвильової ІЧ та теплової ІЧ (див. Рисунок Електромагнітний спектр).

_images/Electromagnetic-Spectrum.png

Електромагнітний спектр

за Victor Blacus (SVG версія File:Electromagnetic-Spectrum.png)

[CC-BY-SA-3.0 (http://creativecommons.org/licenses/by-sa/3.0)]

за посередництвом Wikimedia Commons

http://commons.wikimedia.org/wiki/File%3AElectromagnetic-Spectrum.svg


Взаємодія між сонячною енергією та матеріалами залежить від довжини хвилі; сонячна енергія проходить від Сонця до Землі, а потім до сенсора. Вздовж цього шляха сонячна енергія (NASA, 2013):

  • Проникла - енергія проходить через з зміною у швидкості відповідно до індексу заломлювання для двох середовищ, про які йде мова.

  • Поглинена - енергія передається об’єкту через електронні або молекулярні реакції.

  • Відбита - енергія повертається незміненою з кутом відбивання, що дорівнює куту надходження. Відбивальна здатність це відношення відбитої енергії до тієї, що надійшла до тіла. Довжина хвилі відбитої (а не поглиненої) енергії визначає колір об’єкта.

  • Розсіяна - напрямок поширювання енергії змінюється випадково. Розсіювання Рейлі та Мі - два найбільш важливі типи розсіювання в атмосфері.

  • Випромінена - в дійсності, енергія спочатку поглинається, а потім випромінюється знову, зазвичай на довших довжинах хвиль. Об’єкт розігрівається.

4.1.3. Сенсори

Сенсори можуть знаходитись на борті літака або супутника, вимірюючи електромагнітну радіацію у деяких визначених діапазонах (які зазвичай називаються каналами). В результаті, вимірювання квантуються та перетворюються на цифрове зображення, кожен елемент якого (тобто піксель) має дискретне значення в одиницях цифрових чисел (Digital Number - DN) (NASA, 2013). Результуючі зображення мають різні характеристики (вирізняльні здатності) залежно від сенсора. Виділяють декілька типів вирізняльної здатності:

  • Просторова вирізняльна здатність, зазвичай відповідає роміру пікселя, «це вирізняльна здатність інструмента, необхідна для виокремлення об’єктів, що залежить від розміру детектора, фокусної відстані та висоти сенсора» (NASA, 2013); просторову родільну здатність також називають геометричною роздільною здатністю або IFOV;

  • Спектральна вирізняльна здатність це кількість та положення електромагнітних спектрів (що визначаються двома довжинами хвиль) спектральних каналів (NASA, 2013) багатоспектральних сенсорів, кожному каналу відповідає зображення;

  • Радіометрична вирізняльна здатність, зазвичай вимірюється у бітах (двійкових цифрах), це діапазон можливих значень яскравості, який для знімка відповідає максимальному діапазону DN; наприклад, знімок з вирізняльною здатністю 8 біт має 256 рівнів яскравості (Richards and Jia, 2006);

  • Для супутникових сенсорів, є також часова вирізняльна здатність, яка відповідає часу, необхідному для повторного перегляду тієї самої ділянки Землі (NASA, 2013).

4.1.4. Енергетична світність та відбивальна здатність

Сенсори вимірюють енергетичну світність, яка відповідає яскравості у заданому напрямку до сенсора; також доцільно вирізняти відбивальну здатність як відношення відбитої до загальної енергії.

4.1.5. Спектральна сигнатура

The spectral signature is the reflectance as a function of wavelength (see Figure Криві спектральної відбивальної здатності для чотирьох різних цілей); each material has a unique signature, therefore it can be used for material classification (NASA, 2013).

_images/Spectral_Reflectance_NASA.jpg

Криві спектральної відбивальної здатності для чотирьох різних цілей

(за NASA, 2013)

4.1.6. Земельний покрив

Земельний покрив це матеріал на земній поверхні такий як ґрунт, рослинність, вода, асфальт тощо (Fisher and Unwin, 2005). Кількість та вид класів земельного покриву, що можуть бути ідентифіковані на знімку, можуть істотно різнитися залежно від вирізняльної здатності сенсора, .

4.2. Multispectral satellites

There are several satellites with different characteristics that acquire multispectral images of earth surface. The following satellites are particularly useful for land cover monitoring because images are provided for free and can be downloaded directly from SCP; data have been acquired for the past few decades and the archive is continuously growing with recent images.

4.2.1. Landsat Satellites

Landsat це серія багатоспектральних супутників, розроблених NASA (Національним управлінням з аеронавтики і дослідження космічного простору США), з початку 1970-х.

Знімки Landsat широко вживані у дослідженнях довкілля. Вирізняльні здатності сенсорів Landsat 4 та Landsat 5 наведені в таблиці нижче (за http://landsat.usgs.gov/band_designations_landsat_satellites.php); часова вирізняльна здатність Landsat становить 16 днів (NASA, 2013).

Канали Landsat 4 та Landsat 5

Канали Landsat 4, Landsat 5

Довжина хвилі [мікрометри]

Вирізняльна здатність [метри]

Канал 1 - Синій

0.45 - 0.52

30

Канал 2 - Зелений

0.52 - 0.60

30

Канал 3 - Червоний

0.63 - 0.69

30

Канал 4 - Ближній інфрачервоний (NIR)

0.76 - 0.90

30

Канал 5 - SWIR

1.55 - 1.75

30

Канал 6 - Тепловий інфрачервоний

10.40 - 12.50

120 (передискретизована до 30)

Канал 7 - SWIR

2.08 - 2.35

30

Вирізняльні здатності сенсора Landsat 7 наведено в таблиці нижче (за http://landsat.usgs.gov/band_designations_landsat_satellites.php); часова вирізняльна здатність Landsat становить 16 днів (NASA, 2013).

Канали Landsat 7

Landsat 7 Bands

Довжина хвилі [мікрометри]

Вирізняльна здатність [метри]

Канал 1 - Синій

0.45 - 0.52

30

Канал 2 - Зелений

0.52 - 0.60

30

Канал 3 - Червоний

0.63 - 0.69

30

Канал 4 - Ближній інфрачервоний (NIR)

0.77 - 0.90

30

Канал 5 - SWIR

1.57 - 1.75

30

Канал 6 - Тепловий інфрачервоний

10.40 - 12.50

60 (передискретизована до 30)

Канал 7 - SWIR

2.09 - 2.35

30

Канал 8 - Панхроматичний

0.52 - 0.90

15

Вирізняльні здатності сенсора Landsat 8 наведено в таблиці нижче (за http://landsat.usgs.gov/band_designations_landsat_satellites.php); часова вирізняльна здатність Landsat становить 16 днів (NASA, 2013).

Канали Landsat 8

Landsat 8 Bands

Довжина хвилі [мікрометри]

Вирізняльна здатність [метри]

Канал 1 - Узбережний аерозоль

0.43 - 0.45

30

Канал 2 - Синій

0.45 - 0.51

30

Канал 3 - Зелений

0.53 - 0.59

30

Канал 4 - Червоний

0.64 - 0.67

30

Канал 5 - Ближній інфрачервоний (NIR)

0.85 - 0.88

30

Канал 6 - SWIR 1

1.57 - 1.65

30

Канал 7 - SWIR 2

2.11 - 2.29

30

Канал 8 - Панхроматичний

0.50 - 0.68

15

Канал 9 - Пір’їсті хмари

1.36 - 1.38

30

Канал 10 - Тепловий інфрачервоний (TIRS) 1

10.60 - 11.19

100 (передискретизована до 30)

Канал 11 - Тепловий інфрачервоний (TIRS) 2

11.50 - 12.51

100 (передискретизована до 30)

Великий архів знімків безкоштовно доступний від Геологічної служби США . Для більш докладної інформації щодо безкоштовного завантаження знімків Landsat читайте .

Знімки ідентифікуються траєкторіями та рядками WRS (Глобальна система місцерозташування для Landsat Worldwide Reference System for Landsat ).

4.2.2. Супутник Sentinel-2

Sentinel-2 це багатоспектральний супутник, розроблений Европейським космічним агентством (European Space Agency - ESA) в рамках програми мніторингу земель Copernicus . Sentinel-2 має 13 спектральних каналів з просторовою вирізняльною здатністю 10 м, 20 м та 60 м залежно від каналу, як наведено у таблиці нижче (ESA, 2015).

Канали Sentinel-2

Канали Sentinel-2

Цетральна довжина хвилі [мікрометри]

Вирізняльна здатність [метри]

Канал 1 - Узбережний аерозоль

0.443

60

Канал 2 - Синій

0.490

10

Канал 3 - Зелений

0.560

10

Канал 4 - Червоний

0.665

10

Канал 5 - Червоний край рослинності

0.705

20

Канал 6 - Червоний край рослинності

0.740

20

Канал 7 - Червоний край рослинності

0.783

20

Канал 8 - NIR

0.842

10

Канал 8А - Червоний край рослинності

0.865

20

Канал 9 - Водяна пара

0.945

60

Канал 10 - SWIR - Пір’їсті хмари

1.375

60

Канал 11 - SWIR

1.610

20

Канал 12 - SWIR

2.190

20

Sentinel-2 images are freely available from the ESA website https://scihub.copernicus.eu.

4.2.3. Sentinel-3 Satellite

Sentinel-3 is a satellite developed by the European Space Agency (ESA) in the frame of Copernicus land monitoring services. It carries several instruments, in particular the Ocean and Land Colour Instrument (OLCI) is a push-broom imaging spectrometer acquiring 21 bands in the range 0.4-1.02 μm with a swath width of 1,270km and 300m spatial resolution (ESA, 2013). The revisit time is about 2 days.

Sentinel-3 Bands

Sentinel-3 Bands

Цетральна довжина хвилі [мікрометри]

Oa1

0.400

Oa2

0.4125

Oa3

0.4425

Oa4

0.490

Oa5

0.510

Oa6

0.560

Oa7

0.620

Oa8

0.665

Oa9

0.67375

Oa10

0.68125

Oa11

0.70875

Oa12

0.75375

Oa13

0.76125

Oa14

0.764375

Oa15

0.7675

Oa16

0.77875

Oa17

0.865

Oa18

0.885

Oa19

0.900

Oa20

0.940

Oa21

1.020

4.2.4. Супутник ASTER

Супутник ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer - Покращений космічний радіометр теплового випромінювання та відбивання) був запущений у 1999 за співробітництва між Міністерством міжнародної торігівлі та індустрії Японії (Japanese Ministry of International Trade and Industry - MITI) та NASA. ASTER має 14 каналів вирізняльна здатність яких змінюється відповідно до довжини хвилі: 15 м у видимому та ближньому інфрачервоному, 30 м у короткохвильовому інфрачервоному та 90 м у тепловому інфрачервоному (USGS, 2015). Канали ASTER охарактеризовано в таблиці нижче (через вихід з ладу сенсора дані SWIR, отримані після 1 квітня 2008 недоступні ). Додатковий канал 3B (ближній інфрачервоний оберненого назад огляду) забезпечує стерео покриття.

Канали ASTER

Канали ASTER

Довжина хвилі [мікрометри]

Вирізняльна здатність [метри]

Канал 1 - Зелений

0.52 - 0.60

15

Канал 2 - Червоний

0.63 - 0.69

15

Канал 3N - Ближній інфрачервоний (NIR)

0.78 - 0.86

15

Канал 4 - SWIR 1

1.60 - 1.70

30

Канал 5 - SWIR 2

2.145 - 2.185

30

Канал 6 - SWIR 3

2.185 - 2.225

30

Канал 7 - SWIR 4

2.235 - 2.285

30

Канал 8 - SWIR 5

2.295 - 2.365

30

Канал 9 - SWIR 6

2.360 - 2.430

30

Канал 10 - TIR 1

8.125 - 8.475

90

Канал 11 - TIR 2

8.475 - 8.825

90

Канал 12 - TIR 3

8.925 - 9.275

90

Канал 13 - TIR 4

10.25 - 10.95

90

Канал 14 - TIR 5

10.95 - 11.65

90

4.2.5. Продукти MODIS

The MODIS (Moderate Resolution Imaging Spectroradiometer) is an instrument operating on the Terra and Aqua satellites launched by NASA in 1999 and 2002 respectively. Its temporal resolutions allows for viewing the entire Earth surface every one to two days, with a swath width of 2,330km. Its sensors measure 36 spectral bands at three spatial resolutions: 250m, 500m, and 1,000m (see https://lpdaac.usgs.gov/dataset_discovery/modis).

Доступні декілька продуктів, таких як відбивальність поверхні та вегетаційні індекси. В цьому посібнику ми розглядаємо канали відбивальності поверхні доступні з просторовою вирізняльною здатністю 250 м та 500 м (Vermote, Roger, & Ray, 2015).

Канали MODIS

Канали MODIS

Довжина хвилі [мікрометри]

Вирізняльна здатність [метри]

Канал 1 - Червоний

0.62 - 0.67

250 - 500

Канал 2 - Ближній інфрачервоний (NIR)

0.841 - 0.876

250 - 500

Канал 3 - Синій

0.459 - 0.479

500

Канал 4 - Зелений

0.545 - 0.565

500

Канал 5 - SWIR 1

1.230 - 1.250

500

Канал 6 - SWIR 2

1.628 - 1.652

500

Канал 7 - SWIR 3

2.105 - 2.155

500

Наступні продукти (Версія 6, див. https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table) доступні для завантаження (Vermote, Roger, & Ray, 2015):

  • MOD09GQ: щоденна відбивальність з просторовою вирізняльною здатністю 250 м за Terra MODIS;

  • MYD09GQ: щоденна відбивальність з просторовою вирізняльною здатністю 250 м за Aqua MODIS;

  • MOD09GA: щоденна відбивальність з просторовою вирізняльною здатністю 500 м за Terra MODIS;

  • MYD09GA: щоденна відбивальність з просторовою вирізняльною здатністю 500 м за Aqua MODIS;

  • MOD09Q1: відбивальність з просторовою вирізняльною здатністю 250 м, як композит MOD09GQ (кожний піксель містить найкраще можливе спостереження за 8-денний період);

  • MYD09Q1: відбивальність з просторовою вирізняльною здатністю 250 м, як композит MYD09GQ (кожний піксель містить найкраще можливе спостереження за 8-денний період);

  • MOD09A1: відбивальність з просторовою вирізняльною здатністю 250 м, як композит MOD09GA (кожний піксель містить найкраще можливе спостереження за 8-денний період);

  • MYD09A1: відбивальність з просторовою вирізняльною здатністю 250 м, як композит MYD09GA (кожний піксель містить найкраще можливе спостереження за 8-денний період);

4.2.6. GOES Products

The Geostationary Operational Environmental Satellite-R Series (GOES-R) are geostationary satellites developed for weather monitoring by the National Oceanic and Atmospheric Administration (NOAA) and the NASA (NOAA, 2020).

GOES constellation is composed of GOES-R satellite (also known as GOES-16 that replaced GOES-13 on December 2017), and GOES-S satellite (also known as GOES-17, operational since February 2019). Other satellites (GOES-T and GOES-U) are planned to be launched in the future. For more information please visit https://www.goes-r.gov/mission/mission.html .

GOES geostationary satellites monitor continuously the same area, a very large portion of Earth surface with three geographic coverage regions: Full Disk, Continental United States (CONUS), and Mesoscale. In particular, Full Disk products have hemispheric coverage of 83° local zenith angle, and images are acquired every 5-15 minutes. GOES-16 monitors from 75.2 degrees west longitude, including America, the Atlantic Ocean, and the west coast of Africa. GOES-17 monitors from 137.2 degrees west longitude, including the Pacific Ocean.

GOES sensors include several spectral bands.

GOES Bands

GOES Bands

Цетральна довжина хвилі [мікрометри]

Вирізняльна здатність [метри]

Канал 1 - Синій

0.47

1000

Канал 2 - Червоний

0.64

500

Band 3 - Near Infrared (NIR)

0.87

1000

Band 4 - SWIR - Cirrus

1.38

2000

Канал 5 - SWIR

1.61

1000

Band 6 - SWIR

2.25

2000

4.3. SAR satellites

Synthetic Aperture Radar (SAR) is a technique of active remote sensing that is the sensor platform emits microwaves in order to acquire images of the ground (Richards and Jia, 2006). In fact, the sensor platform emits the radiation (at a specific wavelength) and measures the magnitude and the phase of radiation that bounces back from the ground to the sensor.

Unlike passive sensors, SAR systems can work day and night and can penetrate clouds allowing for the monitoring of surface also with adverse meteorological conditions; depending on the microwave wavelength, the radiation can penetrate different types of materials allowing for different applications (NASA, 2020).

The main SAR systems can be divided according to the wavelength as illustrated in the following table (NASA, 2020):

Main SAR Bands

Канал

Wavelength [centimeters]

Application

X

3.8 – 2.4

High Resolution SAR, urban monitoring, ice and snow, little penetration into vegetation cover

C

7.5 – 3.8

global mapping, change detection, ice, low penetration into vegetation cover

S

15 – 7.5

global mapping, agriculture monitoring, medium penetration into vegetation cover

L

30 – 15

Medium resolution SAR, biomass and vegetation mapping, high penetration into vegetation cover

Usually, SAR sensors can emit and measure different polarizations (i.e. orientation of the microwaves of the electric field), for instance vertical (i.e. polarization oriented in the vertical direction in antenna coordinates) and horizontal (i.e. polarization oriented in the horizontal direction in antenna coordinates) (ESA, 2020).

SAR systems can acquire in both ascending and descending orbits, however the acquired images are affected by the different acquisition geometries, which should be considered when mixing ascending and descending images.

Acquisitions are called swaths and usually they are composed of sub-swaths. With particular acquisition modes, the resolution of pixels along track (the side parallel to the flight direction) can be different than slant-range (the side perpendicular to the flight direction).

SAR phase information is used to perform interferometry (also InSAR) to measure the distance from the sensor to the target (NASA, 2020).

For more information, please read the ESA introduction to SAR and the NASA definition of SAR .

4.3.1. Sentinel-1 Satellites

Sentinel-1 is a Copernicus mission of satellites that operate at C-band to provide SAR imagery at medium resolution (about 10m).

The Sentinel-1 constellation provides high revisit time (about 5 days), a wide swath (250 km), and acquires images in different operational modes. The primary operation mode on land is the Interferometric Wide swath (IW), which is data is acquired in three swaths using the Terrain Observation with Progressive Scanning SAR (TOPSAR) imaging technique (ESA, 2020b).

The Level-1 products systematically delivered by Copernicus are Single Look Complex (SLC, data comprising complex imagery with amplitude and phase) and Ground Range Detected (GRD, Level-1 data with multi-looked intensity only).

Sentinel-1 supports dual polarization, which are horizontal (H) or vertical (V); VV and VH polarimetric channels are available to classify and analyze land cover such as built-up areas or vegetation.

4.4. Класифікація земельного покриву

This chapter provides basic definitions about land cover classifications.

4.4.1. Контрольована класифікація

Напівавтоматична класифікація (також контрольована класифікація) це техніка обробки зображень, яка дозволяє визначити матеріали на знімку відповідно до їх спектральних сигнатур. Існує декілька видів алгоритмів класифікації, але головна мета це створення тематичної карти земельного покриву.

Обробка зобоажень та просторовий ГІС-аналіз потребують спеціалізованого програмного забезпечення такого як Semi-Automatic Classification Plugin для QGIS.

_images/multispectral_classification.jpg

Багатоспектральне зображення оброблене для створення класифікації земельного покриву

(Знімок Landsat надано USGS)

4.4.2. Кольоровий композит

Зазвичай створюється комбінація трьох індивідуальних монохромних зображень, в якій кожному призначається певний колір; така комбінація називається кольоровий композит та корисна для візуальної інтерпретації (NASA, 2013). Кольрові композити можуть бути описані наступним виразом:

«R G B = Br Bg Bb»

де:

  • R відповідає Червоному;

  • G відповідає Зеленому;

  • B відповідає Синьому;

  • Br це номер каналу, що асоціюється з червоним кольором;

  • Bg це номер каналу, що асоціюється з зеленим кольором;

  • Bb це номер каналу, що асоціюється з синім кольором.

Наступний рисунок Кольоровий композит знімка Landsat 8 демонструє кольоровий композит «R G B = 4 3 2» знімка Landsat 8 (для Landsat 7 аналогічний кольровий композит це R G B = 3 2 1; для Sentinel-2 це R G B = 4 3 2) та кольоровий композит «R G B = 5 4 3» (для Landsat 7 аналогічний кольоровий композит це R G B = 4 3 2; для Sentinel-2 це R G B = 8 4 3). Композит «R G B = 5 4 3» корисний для інтерпретації зображень тому що пікселі з рослинністю виглядають червоними (здорова рослинність відбиває значну частину надхідного світла в ближній іфрачервоній зоні, що проявляється у вищих значеннях відбивальності для каналу 5, а відтак і вищих значеннях для пов’язаного червоного кольору).

_images/color_composite.jpg

Кольоровий композит знімка Landsat 8

Дані доступні з Геологічної служби США

4.4.3. Навчальні області

Зазвичай, контрольована класифікація потребує, щоб користувач визначив одну або декілька областей інтересу (Regions of Interest - ROI або навчальних областей) для кожного класу земельного покриву, що визначається на знімку. ROI це полігони окреслені навколо однорідних ділянок зображення, що накладаються на пікселі, які належать до одного класу земельного покриву.

4.4.3.1. Алгоритм нарощування області

Алгоритм нарощування області дозволяє вибрати пікселі подібні до насінини з урахуванням спектральної подібності (тобто спектральної відстані) прилеглих пікселів. В SCP алгоритм нарощування області доступний для створення навчальних областей. Параметр відстань пов’язаний з подібністю значень пікселів (чим нижче значення, тим більш подібні вибрані пікселі) до насінини (вибраної клацанням на пікселі). Додатковий параметр це максимальна ширина, яка є довжиною сторони квадрата з центром в пікселі-насінині, до якого вписана навчальна область (якби всі пікселі мали одні і ті самі значення, то навчальна область також була б квадратною). Мінімальний розмір використовується в якості обмеження (поокремо для кожного каналу) при виборі пікселів, що є більш подібними до насінини, доки їх кількість не досягне принаймні мінімального розміру.

На рисунку Приклад нарощування області центральний пісель використовується в якості насінини (a) для нарощування області одного каналу (b) з параметром спектральної відстані = 0.1; подібні пікселі вибираються для створення навчальної області (c та d).

_images/region_growing.jpg

Приклад нарощування області

4.4.4. Класи та макрокласи

Класи земельного покриву ідентифікуються за довільними кодами ID (тобто унікальними ідентифікаторами). SCP дозволяє призначати ID макрокласу (тобто MC ID) та ID класу (тобто C ID), які є кодами-ідентифікаторами класів земельного покриву. Макроклас це група ROI, що мають різні ID класу, які є зручними за необхідності класифікувати матеріали, що характеризуються різними спектральними сигнатурами, але належать до одного класу земельного покриву. Наприклад, можна віділити траву (ID класу = 1 та ID макрокласу = 1 ) та дерева (ID класу = 2 та ID макрокласу = 1 ) як клас рослинності (ID макрокласу = 1 ). Як показано в таблиці нижче, декілька ID класів можуть бути віднесені до одного і того самого ID макрокласу, але один і той самий ID класу не може бути віднесений до багатьох ID макрокласів.

Приклад макрокласів

Назва макрокласу

ID макрокласу

Назва класу

ID класу

Рослинність

1

Трава

1

Рослинність

1

Дерева

2

Забудова

2

Будівлі

3

Забудова

2

Дороги

4

Відтак, як показано на рис. Приклад макрокласу, класи є підмножинами макрокласу.

_images/macroclass_example.jpg

Приклад макрокласу

Якщо цілі дослідження не вимагають використання макрокласу, тоді один і той самий ID макрокласу може бути призначений всім ROI (наприклад, ID макрокласу = 1) та значення макрокласу не братимуться до уваги в процесі класифікації.

4.4.5. Алгоритми класифікації

Спектральні сигнатури (спектральні характеристики) визначених класів земельного покриву розраховуються з урахуванням значень пікселів кожної ROI, що має той самий ID класу (або ID макрокласу). Відтак, алгоритм класифікації класифікує все зображення шляхом співставлення спектральних характеристик кожного пікселя з спектральними характеристиками визначених класів земельного покриву. SCP здатний реалізовувати наступні алгоритми класифікації.

4.4.5.1. Мінімальної відстані

Алгоритм мінімальної відстані розраховує Евклідову відстань \(d(x, y)\) між спектральними сигнатурами пікселів зображення та навчальними спектральними сигнатурами за наступною формулою:

\[d(x, y) = \sqrt{ \sum_{i=1}^{n} (x_i - y_i)^2}\]

де:

  • \(x\) = вектор спектральної сигнатури пікселя зображення;

  • \(y\) = вектор спектральної сигнатури навчальної області;

  • \(n\) = кількість каналів знімка.

Відтак, відстань розраховується для кожного пікселя на знімку з присвоєнням класу найближчої спектральної сигнатури відповідно до наступної дискримінантної функції (з змінами за Richards and Jia, 2006):

\[x \in C_k \iff d(x, y_k) < d(x, y_j) \forall k \neq j\]

де:

  • \(C_k\) = клас земельного покриву \(k\);

  • \(y_k\) = спектральна сигнатура класу \(k\);

  • \(y_j\) = спектральна сигнатура класу \(j\).

Є можливість визначення порогової величини \(T_i\) з метою виключення з класифікації пікселів, що знаходяться нижче цього значення:

\[ \begin{align}\begin{aligned}x \in C_k \iff d(x, y_k) < d(x, y_j) \forall k \neq j\\and\\d(x, y_k) < T_i\end{aligned}\end{align} \]

4.4.5.2. Максимальної вірогідності

Maximum Likelihood algorithm calculates the probability distributions for the classes, related to Bayes’ theorem, estimating if a pixel belongs to a land cover class. In particular, the probability distributions for the classes are assumed the of form of multivariate normal models (Richards & Jia, 2006). In order to use this algorithm, a sufficient number of pixels is required for each training area allowing for the calculation of the covariance matrix. The discriminant function, described by Richards and Jia (2006), is calculated for every pixel as:

\[g_k(x) = \ln p(C_k) - \frac{1}{2} \ln | \Sigma_{k} | - \frac{1}{2} (x - y_k)^t \Sigma_{k}^{-1} (x - y_k)\]

де:

  • \(C_k\) = клас земельного покриву \(k\);

  • \(x\) = вектор спектральної сигнатури пікселя зображення;

  • \(p(C_k)\) = імовірність, що вірний клас це \(C_k\);

  • \(| \Sigma_{k} |\) = визначник матриці коваріації даних в класі \(C_k\);

  • \(\Sigma_{k}^{-1}\) = обернена матриця коваріації;

  • \(y_k\) = вектор спектральної сигнатури класу \(k\).

Таким чином:

\[x \in C_k \iff g_k(x) > g_j(x) \forall k \neq j\]
_images/maximum_likelihood.jpg

Приклад максимальної вірогідності

Також є можливість визначення порогової величини дискрмінантної функції з метою виключення з класифікації пікселів, що знаходяться нижче цього значення. З урахуванням порогової величини \(T_i\) умова класифікації набуває вигляду:

\[ \begin{align}\begin{aligned}x \in C_k \iff g_k(x) > g_j(x) \forall k \neq j\\and\\g_k(x) > T_i\end{aligned}\end{align} \]

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

4.4.5.3. Картографування спектрального кута

Алгоритм картографування спектрального кута розраховує спектральний кут між спектральними сигнатурами пікселів зображення та навчальними спектральними сигнатурами. Спектральний кут \(\theta\) визначається як (Kruse et al., 1993):

\[\theta(x, y) = \cos^{-1} \left( \frac{ \sum_{i=1}^{n} x_i y_i } { \left( \sum_{i=1}^{n} x_i^2 \right)^\frac{1}{2} * \left( \sum_{i=1}^{n} y_i^2 \right)^\frac{1}{2} } \right)\]

Де:

  • \(x\) = вектор спектральної сигнатури пікселя зображення;

  • \(y\) = вектор спектральної сигнатури навчальної області;

  • \(n\) = кількість каналів знімка.

Відтак піксель належить до класу, що має найменший кут, тобто:

\[x \in C_k \iff \theta(x, y_k) < \theta(x, y_j) \forall k \neq j\]

де:

  • \(C_k\) = клас земельного покриву \(k\);

  • \(y_k\) = спектральна сигнатура класу \(k\);

  • \(y_j\) = спектральна сигнатура класу \(j\).

_images/spectral_angle_mapping.jpg

Приклад картографування спектрального кута

З метою виключення з класифікації пікселів нижче цього значення можливо призначити порогову величину \(T_i\):

\[ \begin{align}\begin{aligned}x \in C_k \iff \theta(x, y_k) < \theta(x, y_j) \forall k \neq j\\and\\\theta(x, y_k) < T_i\end{aligned}\end{align} \]

Алгоритм картографування спектрального кута широко застосовується, особливо з гіперспектральними даними.

4.4.5.4. Класифікація паралелепіпеда

Класифікація паралелепіпеда це алгоритм, який враховує діапазон значень для кожного каналу, формуючи багатовимірний паралелепіпед, що визначає клас земельного покриву. Піксель відноситься до класу, якщо його значення знаходяться всередині паралелепіпеда. Одним з основних обмежень є те, що пікселі, сигнатури яких знаходяться в областях перекриття двох або більше паралелепіпедів, не можуть бути класифіковані (Richards and Jia, 2006).

4.4.5.5. Класифікація сигнатур земельного покриву

This classification allows for the definition of spectral thresholds for each training input signature (a minimum value and a maximum value for each band). The thresholds of each training input signature define a spectral region belonging to a certain land cover class.

Спектральні сигнатури пікселів зображення порівнюються з спектральними навчальними сигнатурами; піксель належить до класу X, якщо його спектральна сигнатура повністю міститься в спектральній області, що визначається класом X. У випадку, якщо пікселі потрапляють до області перекриття або поза межі будь-якої спектральної області, є можливість застосувати додаткові алгоритми класифікації (наприклад, Мінімальної відстані, Максимальної вірогідності, Картографування спектрального кута) з урахуванням спектральних характеристик первинних входових сигнатур.

Рисунок нижче схематично ілюструє Класифікація сигнатур земельного покриву для простого випадку двох спектральних каналів \(x\) та \(y\). Визначені користувачем спектральні області означують три класи (\(g_a\), \(g_b\) та \(g_c\)). Точка \(p_1\) належить до класу \(g_a\), а точка \(p_2\) - до класу \(g_b\). Однак точка \(p_3\) потрапляє всередину спектральних областей обох класів \(g_b\) та \(g_c\) (області, що перекриваються); в такому випадку тока \(p_3\) залишиться некласифікованою або буде класифікована за додатковим алгоритмом класифікації. Точка \(p_4\) знаходиться поза межами будь-якої спектральної області, тому вона залишиться некласифікованою або буде класифікована за додатковим алгоритмом класифікації. За умови, що точка \(p_4\) належить до класу \(g_c\), спектральна область може бути розширена, щоб включити точку \(p_4\) .

_images/land_cover_signature_classification.jpg

Класифікація сигнатур земельного покриву

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

На рисунку Графік спектральних діапазонів показано спектральні діапазони трьох класів (\(g_a\), \(g_b\) та \(g_c\)); кольорові лінії всередині діапазонів (тобто напівпрозорої області) представляють спектральні сигнатури пікселів, що визначають верхній та нижній ліміти відповідних діапазонів. Піксель \(p_1\) (крапкова лінія) належить до класу \(g_b\) тому що його спектральна сигнатура знаходиться повністю всередині діапазону класу \(g_b\) (у верхньому ліміті); піксель \(p_2\) (пунктирна лінія) не класифікований, тому що його спектральна сигнатура не входить повністю до жодного діапазону; піксель \(p_3\) (крапкова лінія) належить до класу \(g_a\).

_images/land_cover_signature_classification_plot.jpg

Графік спектральних діапазонів

Варто зазначити, що ці спектральні порогові величини можуть бути застосовані до будь-якої сигнатури, безвідносно її спектральних характеристик; ця функція може бут дуже корисною для відокремлення подібних спектральних сигнатур, що відрізняються лише в одному каналі, з визначенням порогових величин, які включають або виключають конкретні сигнатури. Фактично, класи відокремлюються коректно якщо їх спектральні області не перекриваються принаймні в одному каналі. Звичайно, навіть за наявності перекриття спектральних областей є можливість, що жоден піксель не потрапить до області перекриття і не буде невірно класифікований; верхній (або нижній) ліміт області не передбачає існування на зображенні будь-якої спектральної сигнатури, що має максимальне (або мінімальне) значення діапазону за всіма каналами (наприклад піксель math:p_1 рисунка Графік спектральних діапазонів не міг би існувати).

Однією з головних переваг Класифікації сигнатур земельних покривів є можливість вибору пікселів та включення їх сигнатур до спектрального діапазону; відтак, класифікація повинна бути безпосереднім представленням класу, що очікується для кожної спектральної сигнатури. Це дуже доречно для класифікації єдиного класу земельного покриву (визначеного за специфічними спектральними пороговими величинами) та залишає некласифікованою іншу частину зображення, яка не представляє інтересу для цілей класифікації.

4.4.5.6. Растр алгоритму

Растр алгоритму представляє «відстань» (відповідно до визначення алгоритму класифікації) пікселя зображення до певної спектральної сигнатури.

Здебільшого растр алгоритму створюється для кожної спектральної сигнатури, використаної в якості навчальних входових даних. Значення кожного пікселя це результат розрахунків алгоритму для конкретної спектральної сигнатури. Відтак, піксель належить до класу X якщо значення растра алгоритму, що відповідає класу X є найнижчим у випадку Мінімальної відстані або Картографування спектрального кута (або найвищим у випадку Максимальної вірогідності).

Залежно від класифікації може бути зроблена комбінація растрів алгоритму, щоб створити растр найменших «відстаней» (тобто пікселі мають значення растра алгоритму, що співвідноситься з класом до якого вони належать за класифікацією). Таким чином, цей растр може бути корисним для визначення пікселів, що потребують додаткового збору подібних спектральних сигнатур (див. Classification preview).

4.4.6. Machine Learning

Machine Learning is a broad set of classification techniques that aim to build mathematical models based on training data.

In general, Machine Learning algorithms split the data in (ESA, 2019):

  • Training Dataset: the sample data used to fit the model;

  • Validation Dataset: the sample data used to tune the model parameters to fit on the training dataset;

  • Test Dataset: the sample data used to provide an evaluation of the final model;

Usually, the training and model evaluation are performed iteratively.

4.4.6.1. Random Forest

Random Forest is a particular machine learning technique, based on the iterative and random creation of decision trees (i.e. a set of rules and conditions that define a class).

First, the input features should be defined, which can be spectral bands or ancillary rasters. Навчальні області should be created to define the classes used for training the model.

Random Forest calculates several random decision trees, based on the following parameters:

  • number of training samples: is the number of training data (pixels) randomly used to train the model; it should be set lower than total training input pixels;

  • number of trees: is the number of decision trees; the more the number of trees, the more is the model accuracy, but it also increases the calculation time.

For instance, a decision tree could be defined as:

  • class 1 = band 1 > 0.1 –> band 2 < 0.3 –> band 3 > 0.4

  • class 2 = band 1 > 0.4 –> band 2 > 0.6 –> band 3 < 0.1

  • class 3 = band 1 < 0.7 –> band 2 > 0.1 –> band 3 < 0.5

Random Forest creates several decision trees randomly. Usually, the Gini coefficient is calculated to split the trees. Therefore, a model based on the decision trees is created and used to classify all the pixels.

A pixel is classified according to the majority vote of decision trees, for example a pixel is classified as class 1 if most decision trees evaluated it as class 1. Also, a confidence layer is produced, which measures the uncertainty of the model based on training data.

Random Forest can be used to evaluate the importance of input features, according to the contribution thereof to the model.

4.4.6.2. Multi-Layer Perceptron

Multi-Layer Perceptron is a supervised classification algorithm that is based on the definition of an artificial neural network. A Multi-Layer Perceptron is made of an input layer, one or more hidden layers (made of a defined number of neurons that are fully connected by non-linear activation functions), and the output layer (also read this documentation

Several parameters can be defined as described here

4.4.6.3. Support Vector Machine

Support Vector Machine is a supervised classification algorithm that is based on the calculation of hyperplanes in order to separate input data values.

Several parameters can be defined as described at this link

4.4.7. Спектральна відстань

Для того, щоб оцінити, чи подібні класи настільки, що це може призвести до помилок класифікації, доцільно оцінювати спектральну відстань (або відокремність) між навчальними сигнатурами та пікселями. В SCP реалізовано наступні алгоритми для оцінки подібності спектральних сигнатур.

4.4.7.1. Відстань Джефріса-Мацусіти

Відстань Джефріса-Мацусіти розраховує відокремність двох розподілів імовірностей. Це може бути особливо змістовно для оцінювання результатів класифікацій Максимальної вірогідності.

Відстань Джефріса-Мацусіти \(J_{xy}\) розраховується як (Richards and Jia, 2006):

\[J_{xy} = 2 \left( 1 - e^{-B} \right)\]

де:

\[B = \frac{1}{8} (x - y)^t \left( \frac{\Sigma_{x} + \Sigma_{y}}{2} \right)^{-1} (x - y) + \frac{1}{2} \ln \left( \frac{ | \frac{ \Sigma_{x} + \Sigma_{y}}{2} | }{ | \Sigma_{x} |^{\frac{1}{2}} | \Sigma_{y} |^{\frac{1}{2}} } \right)\]

де:

  • \(x\) = вектор першої спектральної сигнатури;

  • \(y\) = вектор другої спектральної сигнатури;

  • \(\Sigma_{x}\) = матриця коваріації вибірки \(x\);

  • \(\Sigma_{y}\) = матриця коваріації вибірки \(y\);

Відстань Джефріса-Мацусіти асимптотна до 2, коли сигнатури абсолютно різні, і наближається до 0, коли сигнатури ідентичні.

4.4.7.2. Спектральний кут

Спектральний кут найбільш придатний для оцінки алгоритму Картографування спектрального кута. Спектральний кут :theta визначається як (Kruse et al., 1993):

\[\theta(x, y) = \cos^{-1} \left( \frac{ \sum_{i=1}^{n} x_i y_i } { \left( \sum_{i=1}^{n} x_i^2 \right)^\frac{1}{2} * \left( \sum_{i=1}^{n} y_i^2 \right)^\frac{1}{2} } \right)\]

Де:

  • \(x\) = вектор спектральної сигнатури пікселя зображення;

  • \(y\) = вектор спектральної сигнатури навчальної області;

  • \(n\) = кількість каналів знімка.

Спектральний кут іде від 0, коли сигнатури ідентичні, до 90, коли сигнатури абсолютно різні.

4.4.7.3. Евклідова відстань

Евклідова відстань особливо корисна для оцінювання результатів класифікацій Мінімальної відстані. Фактично відстань визначається як:

\[d(x, y) = \sqrt{ \sum_{i=1}^{n} (x_i - y_i)^2}\]

де:

  • \(x\) = вектор першої спектральної сигнатури;

  • \(y\) = вектор другої спектральної сигнатури;

  • \(n\) = кількість каналів знімка.

Евклідова відстань дорівнює 0 якщо сигнатури ідентичні і зростає із зростанням спектральної відстані між сигнатурами.

4.4.7.4. Подібність Брея-Кертіса

Подібність Брея-Кертіса це статистика, що використовується для оцінювання спорідненості між двома вибірками (читайте). Вона корисна загалом для оцінки подібності спектральних сигнатур, подібність Брея-Кертіса \(S(x, y)\) розраховується як:

\[S(x, y) = 100 - \left( \frac{\sum_{i=1}^{n} | (x_i - y_i) |}{\sum_{i=1}^{n} x_i + \sum_{i=1}^{n} y_i} \right) * 100\]

де:

  • \(x\) = вектор першої спектральної сигнатури;

  • \(y\) = вектор другої спектральної сигнатури;

  • \(n\) = кількість каналів знімка.

Подібність Брея-Кертіса розраховується у відсотках і змінюється від 0, коли сигнатури абсолютно різні, до 100, коли спектральні сигнатури ідентичні.

4.4.8. Результат класифікації

Результатом процесу класифікації є растр (дивіться приклад класифікації Landsat на рисунку Класифікація Landsat), де значення пікселів співвідносяться з унікальними ідентифікаторами класів і кожний колір представляє клас земельного покриву.

_images/Landsat_classification.jpg

Класифікація Landsat

Дані доступні з Геологічної служби США

Певна кількість помилок може таплятись в класифікації земельного покриву (тобто пікселі віднесені до помилкового класу земельного покриву) через спектральну подібність класів або невірне визначення класу під час збору ROI.

4.4.9. Оцінка точності

Після процесу класифікації доцільно оцінити точність класифікації земельного покриву з метою ідентифікації та вимірювання похибок карти. Зазвичай оцінка точності виконується з розрахунком матриці похибок, яка є таблицею, що порівнює інформацію карти з контрольними даними (тобто з даними підсупутникових спостережень) для певної кількості пробних ділянок (Congalton and Green, 2009)..

Наступна таблиця схематично зображує матрицю похибок, де k це число класів, що ідентифікуються в класифікації земельного покриву, а n це загальна кількість зібраних елементів вибірки. Елементи основної діагоналі (aii) це кількість вірно ідентифікованих елементів, а інші елементи являють собою похибку класифікації.

Схема матриці похибок

Підсупутникове спостереження 1

Підсупутникове спостереження 2

Підсупутникове спостереження k

Сума

Клас 1

\(a_{11}\)

\(a_{12}\)

\(a_{1k}\)

\(a_{1+}\)

Клас 2

\(a_{21}\)

\(a_{22}\)

\(a_{2k}\)

\(a_{2+}\)

Клас k

\(a_{k1}\)

\(a_{k2}\)

\(a_{kk}\)

\(a_{k+}\)

Сума

\(a_{+1}\)

\(a_{+2}\)

\(a_{+k}\)

\(n\)

Відтак, можливо розрахувати загальну точність як відношення кількості елементів, що класифіковані вірно (сума основної діагоналі), до загальної кількості елементів вибірки n (Congalton and Green, 2009).

The overall accuracy (also expressed in percentage) is defined as:

\[O = \sum_{i=1}^{k} a_{ii} / n\]

The user’s accuracy for each class is defined as the ratio (also expressed in percentage) between correct samples and the row total:

\[U_i = a_{ii} / a_{i+}\]

The commission error \(CE_i = 1 - U_i\) corresponds to pixels classified as class \(i\) that actually belong to a different class.

The producer’s accuracy for each class is calculated as the ratio (also expressed in percentage) between correct samples and the column total:

\[P_i = a_{ii} / a_{+i}\]

The omission error \(OE_i = 1 - P_i\) corresponds to pixels actually belonging to class \(i\) that were classified erroneously as a different class.

It is recommended to calculate the area based error matrix (Olofsson et al., 2014) where each element represents the estimated area proportion of each class. This allows for estimating the unbiased user’s accuracy and producer’s accuracy, the unbiased area of classes according to reference data, and the standard error of area estimates.

For further information, the following documentation is freely available: Landsat Data Users Handbook.

4.5. Image processing

Remote sensing images can be processed in various ways in order to obtain classification, indices, or other derived information that can be useful for land cover characterization.

4.5.1. Аналіз головних компонент

Аналіз головних компонент (Principal Component Analysis - PCA) це метод зменшення вимірності змінних (каналів) до головних компонент (JARS, 1993).

Трансформація головних компонент надає новий набір каналів (головних компонент), які мають наступні характеристики: головні компоненти не корелюють; кожна наступна компонента має дисперсію меншу, ніж попередня компонента. Відтак, це ефективний метод виокремлення інформації та ущільнювання даних (Ready and Wintz, 1973).

У випадку знімка з N спектральними каналами головні компоненти отримуються розрахунком матриці (Ready and Wintz, 1973; Richards and Jia, 2006):

\[Y = D^t X\]

де:

  • \(Y\) = вектор головних компонент

  • \(D\) = матриця власних векторів матриці коваріації \(C_x\) в просторі X

  • \(t\) позначає транспонування вектора

А \(X\) розраховується як:

\[X = P - M\]
  • \(P\) = вектор спектральних значень, що відповідають кожному пікселю

  • \(M\) = вектор середніх значень, що відповідають кожному каналу

Відтак, середнє \(X\), що відповідає кожному каналу становить 0. \(D\) формується власними векторами (матриці коваріації \(C_x\)), впорядкованими як власні числа від максимуму до мініму, для отримання максимальної дисперсії у першій компоненті. Таким чином, головні компоненти не корелюють та кожна наступна компонента має дисперсію меншу, ніж попередня (Ready and Wintz, 1973).

Зазвичай перші дві компоненти містять понад 90% дисперсії. Наприклад, перша головна компонента може бути відображена у Кольоровий композит для підкреслення класів Земельний покрив або використана в якості входових даних для Контрольована класифікація.

4.5.2. Панхроматичне об’єднання

Панхроматичне об’єднання або пан-шарпенинг це поєднання спектральної інформації багатоспектральних каналів (MS), які мають нижчу вирізняльну здатність (для каналів Landsat просторова вирізняльна здатність становить 30 м), з просторовою вирізняльною здатністю панхроматичного каналу (PAN), яка для Landsat 7 та 8 становить 15 м. Результатом є багатоспектральне зображення з просторовою вирізняльною здатністю панхроматичного каналу (тобто 15 м). В SCP застосовується перетворення Бровея, де перетворені значення кожного багатоспектрального каналу розраховуються як (Johnson, Tateishi and Hoan, 2012):

\[MSpan = MS * PAN / I\]

де \(I\) це Інтенсивність, яка є функцією багатоспектральних каналів.

Відповідно до декількох тестів, проведених з застосуванням SCP, для I визначено наступні вагові коефіцієнти. Для Landsat 8, Інтенсивність розраховується як:

\[I = (0.42 * Blue + 0.98 * Green + 0.6 * Red ) / 2\]

Для Landsat 7, Інтенсивність розраховується як:

\[I = (0.42 * Blue + 0.98 * Green + 0.6 * Red + NIR) / 3\]
_images/pan_sharpening_comparison.jpg

Приклад знімка Landsat 8, який був підданим панхроматичному об’єднанню. Ліворуч первинні багатоспектральні канали (30 м); праворуч панхроматично об’єднані (15 м)

Дані доступні з Геологічної служби США

4.5.3. Спектральні індекси

Спектральні індекси це математичні дії між спектральними каналами спрямовані на отримання інформації про рослинний покрив (JARS, 1993). Один з найбільш популярних індексів це вегетаційний індекс нормалізованої різниці (англ. Normalized Difference Vegetation Index - NDVI), що визначається як (JARS, 1993):

\[NDVI = ( NIR - Red ) / ( NIR + Red )\]

Значення NDVI варіюють від -1 до 1. Густа та здорова рослинність демонструє вищі значення, а ділянки без рослинного покриву характеризуються низькими значеннями NDVI.

Інший індекс це підсилений вегетаційний індекс (англ. Enhanced Vegetation Index - EVI), який намагається враховувати атмосферні ефекти, такі як енергетична світність, відбита від атмосфери, розраховуючи різницю між синім та червоним каналами (Didan,et al., 2015). EVI визначається як:

\[EVI = G ( NIR - Red ) / ( NIR + C_1 Red - C_2 Blue + L)\]

де: \(G\) масштабний коефіцієнт, \(C_1\) та \(C_2\) коефіцієнти атмосферних ефектів та \(L\) фактор для урахування диференційного NIR та червоного випромінного передавання через рослинний полог. Типові значення коефіцієнтів становлять: \(G = 2.5\), \(L = 1\), \(C_1 = 6\), \(C_2 = 7.5\) (Didan,et al., 2015).

4.5.4. Clustering

Clustering is the grouping of pixels based on spectral similarity (e.g. Евклідова відстань or Спектральний кут) calculated for a multispectral image (Richards and Jia, 2006).

Clustering can be used for unsupervised classification or for the automatic selection of spectral signatures. It is worth noticing that, while Контрольована класифікація produces a classification whit the classes identified during the training process, the classes produced by clustering (i.e. clusters) have no definition and consequently the user must assign a land cover label to each class.

The main advantage of clustering resides in automation. Of course, clusters do not necessarily represent a particular land cover type and additional processing could be required for producing an accurate classification.

There are several types of clustering, mainly based on iterative methods; the following are the algorithms provided in SCP.

4.5.4.1. K-means

The K-means method is based on the calculation of the average spectral signature of clusters (Wikipedia, 2017; JARS, 1993).

At first, the user defines the number of clusters expected in the image, which correspond to as many spectral signatures (i.e. seeds). Starting spectral signatures can be selected in various ways (e.g. randomly, provided by the user, calculated automatically from image values).

During the first iteration clusters are produced calculating the pixel spectral distance with initial spectral signatures. The algorithms Евклідова відстань or Спектральний кут can be used for distance calculation. Pixels are assigned according to the most similar spectral signature, therefore producing clusters.

Then, the average spectral signature is calculated for each cluster of pixels, resulting in the spectral signatures that will be used in the following iteration.

This process continues iteratively producing clusters and mean spectral signatures, until one of the following condition is verified:

  • the spectral distance between the spectral signatures produced in this iteration with the corresponding ones produced in the previous iteration is lower than a certain threshold;

  • the maximum number of iterations is reached.

After the last iteration, a raster of clusters is produced using the spectral signatures derived from the last iteration.

4.5.4.2. ISODATA

The ISODATA (Iterative Self-Organizing Data Analysis Technique) method is similar to K-means but with the additional steps of merging clusters having similar spectral signatures and splitting clusters having too high variability (i.e. standard deviation) of spectral signatures (Ball & Hall, 1965). Following, the SCP implementation of ISODATA is described.

At first, the user defines the number of clusters expected in the image, which correspond to as many spectral signatures (i.e. seeds). Starting spectral signatures can be selected in various ways (e.g. randomly, provided by the user, calculated automatically from image values). Initial parameters provided by user are:

  • \(C\) = number of desired clusters

  • \(N_{min}\) = minimum number of pixels for a cluster

  • \(\sigma_{t}\) = maximum standard deviation threshold for splitting

  • \(D_{t}\) = distance threshold for merging

During the first iteration clusters are produced calculating the Евклідова відстань of pixels with initial spectral signatures. Pixels are assigned according to the most similar spectral signature, therefore producing clusters.

Therefore, the following parameters are calculated:

  • \(N_{i}\) = number of pixels of cluster \(i\)

  • \(S_{i}\) = average spectral signature of cluster \(i\)

  • \(AVERAGEDIST_{i}\) = average distance of cluster \(i\) with the seed spectral signature

  • \(AVERAGEDISTANCE\) = overall average distance of all clusters

  • \(\sigma_{ij}\) = standard deviation of cluster \(i\) in band \(j\)

  • \(\sigma max_{i}\) = maximum standard deviation of cluster \(i\) (i.e. \(max( \sigma_{ij} )\))

  • \(k_{i}\) = band where \(\sigma max_{i}\) occurred

  • \(Sk_{i}\) = value of \(S_{i}\) at band \(k_{i}\)

  • \(P\) = number of clusters

Then, for each cluster \(i\), if \(N_{i}\) < \(N_{min}\) , then the cluster \(i\) is discarded.

If \(P\) <= \(C\) then try to split clusters. For each cluster \(i\):

  • If \(\sigma max_{i}\) > \(\sigma_{t}\) :

    • If ((\(AVERAGEDIST_{i}\) > \(AVERAGEDISTANCE\)) AND (\(N_{i}\) > (2 * \(N_{min}\) + 2) )) OR (\(C\) > 2 * \(P\)):

      • create a new spectral signature \(S_{p + 1}\) = \(S_{i}\)

      • in \(S_{i}\) set the value \(Sk_{i}\) = \(Sk_{i}\) + \(\sigma max_{i}\)

      • in \(S_{p + 1}\) set the value \(Sk_{p + 1}\) = \(Sk_{i}\) - \(\sigma max_{i}\)

      • \(P\) = \(P\) + 1

      • start a new iteration

If \(P\) > (2 * \(C\)) then try to merge clusters.

  • For each combination \(xy\) of spectral signatures calculate \(D_{xy}\) = Евклідова відстань of spectral signatures \(S_{x}\) and \(S_{y}\) .

  • If the minimum \(D_{xy}\) is greater than \(D_{t}\):

    • S_{i} = (\(N_{i}\) * S_{i} + \(N_{j}\) * S_{j})/(\(N_{i}\) + \(N_{j}\))

    • discard S_{j}

    • \(P\) = \(P\) - 1

    • start a new iteration

After the last iteration, a raster of clusters is produced using the spectral signatures derived from the last iteration. The number of clusters can vary according to the processes of splitting and merging.


4.6. Перерахунок знімка у значення відбивальності

В цьому розділі наведено інформацію щодо способу перерахунку у значення відбивальності, реалізованому у SCP.

4.6.1. Енергетична світність на апертурі сенсора

Енергетична світність це “потік енергії (переважно випромінної або надхідної) на одиницю просторового кута поверхні, що залишає одиницю площі поверхні в заданому напрямку”, “Енергетична світність вимірюється сенсором та певною мірою залежить від відбивальності” (NASA, 2011, p. 47).

Знімки, такі як Landsat або Sentinel-2, складаються з декількох каналів та файлу метаданих, який містить інформацію необхідну для перерахунку у значення відбивальності.

Landsat images are provided in radiance, scaled prior to output. For Landsat images Spectral Radiance at the sensor’s aperture (\(L_{\lambda}\), measured in [watts/(meter squared * ster * \(\mu m\))]) is given by (https://www.usgs.gov/core-science-systems/nli/landsat/using-usgs-landsat-level-1-data-product):

\[L_{\lambda} = M_{L} * Q_{cal} + A_{L}\]

де:

  • \(M_{L}\) = залежний від каналу множильний коефіцієнт перемасштабовування з метаданих Landsat (RADIANCE_MULT_BAND_x, де x це номер каналу)

  • \(A_{L}\) = залежний від каналу адитивний коефіцієнт перемасштабовування з метаданих Landsat (RADIANCE_ADD_BAND_x, де x це номер каналу)

  • \(Q_{cal}\) = дискретизовані та калібровані значення пікселів стандартного продукту (DN)

Знімки Sentinel-2 (Level-1C) постачаються вже попередньо змасштабованими у Відбивальність на поверхні атмосфери (ТОА) (ESA, 2015).

4.6.2. Відбивальність на поверхні атмосфери (ТОА)

Знімки у значеннях енергетичної світності можуть бути перераховані у відбивальність на поверхні атмосфери (ТОА) (комбіновану відбивальність земної поверхні та атмосфери) з метою зменшення мінливості між сценами шляхом нормування значень енергетичної освітленості сонцевим проміннямю. Відбивальність ТОА (\(\rho_{p}\)), яка є безрозмірним відношенням відбитої до загальної потужності енергії (NASA, 2011), розраховується як:

\[\rho_{p} = (\pi * L_{\lambda} * d^{2} )/ (ESUN_{\lambda} * cos\theta_{s})\]

де:

  • \(L_{\lambda}\) = спектральна густина енергетичної світності на апертурі сенсора (енергетична світність на супутнику)

  • \(d\) = Earth-Sun distance in astronomical units (provided with Landsat 8 metadata file, and an excel file is available from http://landsathandbook.gsfc.nasa.gov/excel_docs/d.xls )

  • \(ESUN_{\lambda}\) = середня екзоатмосферна енергетична освітленість сонцевим промінням

  • \(\theta_{s}\) = сонцевий зенітний кут в градусах, який дорівнює \(\theta_{s}\) = 90° - \(\theta_{e}\), де \(\theta_{e}\) це висота Сонця

Варто зауважити, що дані Landsat 8 постачаються з залежними від каналу коефіцієнтами перемасштабовування, які дозволяють здійснювати безпосереднє перерахування з DN до відбивальності TOA.

Sentinel-2 images are already provided in scaled TOA reflectance, which can be converted to TOA reflectance with a simple calculation using the Quantification Value provided in the metadata (see https://sentinel.esa.int/documents/247904/349490/S2_MSI_Product_Specification.pdf ).

Sentinel-3 images are already provided in scaled TOA radiance. Conversion to reflectance is performed applying the coefficients scale_factor and add_offset provided in the metadata of each band. The ancillary raster tie_geometries.nc provides the value of sun zenith angle and the ancillary raster instrument_data provides information about the solar flux for each band, which are used for the conversion to reflectance with the correction for sun angle. In addition, the georeferencing of the bands is performed using the ancillary raster geo_coordinates.nc which provides coordinates of every pixel.

4.6.3. Відбивальність поверхні

Для вимірювання відбивальності на земній поверхні повинен бути врахований вплив атмосфери (тобто збурення відбивальності, яке залежить від довжини хвилі).

Відповідно до Moran et al. (1992), відбивальність земної поверхні (\(\rho\)) це:

\[\rho = [\pi * (L_{\lambda} - L_{p}) * d^{2}]/ [T_{v} * ( (ESUN_{\lambda} * cos\theta_{s} * T_{z} ) + E_{down} )]\]

де:

  • \(L_{p}\) енергетична світність, відбита від атмосфери

  • \(T_{v}\) пропускальна здатність атмосфери в напрямку огляду

  • \(T_{z}\) пропускальна здатність атмосфери в напрямку освітлення

  • \(E_{down}\) низхідна енергетична освітленість розсіяним промінням

Таким чином, для розрахунку \(\rho\) необхідно провести декілька атмосферних вимірювань (з метою отримання поправок на підставі абсолютних фізичних величин). В якості альтернативи можна скористатись відносними техніками на основі знімка, які не передбачають проведення позалабораторних вимірювань під час отримання знімка. Варто зазначити, що для даних Landsat 8 доступні Surface Reflectance High Level Data Products (для більш докладної інформації читайте http://landsat.usgs.gov/CDR_LSR.php).

4.6.4. Корекція DOS1

The Dark Object Subtraction (DOS) is a family of image-based atmospheric corrections. Chavez (1996) explains that «the basic assumption is that within the image some pixels are in complete shadow and their radiances received at the satellite are due to atmospheric scattering (path radiance). This assumption is combined with the fact that very few targets on the Earth’s surface are absolute black, so an assumed one-percent minimum reflectance is better than zero percent”. It is worth pointing out that the accuracy of image-based techniques is generally lower than physically-based corrections, but they are very useful when no atmospheric measurements are available as they can improve the estimation of land surface reflectance. The path radiance is given by (Sobrino et al., 2004):

\[L_{p} = L_{min} - L_{DO1\%}\]

де:

  • \(L_{min}\) = «radiance that corresponds to a digital count value for which the sum of all the pixels with digital counts lower or equal to this value is equal to the 0.01% of all the pixels from the image considered” (Sobrino et al., 2004, p. 437), therefore the radiance obtained with that digital count value (\(DN_{min}\))

  • \(L_{DO1\%}\) = енергетична світність темного об’єкта, що за припущенням має значення відбивальності 0.01

Зокрема для знімків Landsat:

\[L_{min} = M_{L} * DN_{min} + A_{L}\]

Знімки Sentinel-2 перераховуються у значення енергетичної світності до проведення розрахунків DOS1.

The radiance of Dark Object is given by (Sobrino et al., 2004):

\[L_{DO1\%} = 0.01 * [(ESUN_{\lambda} * cos\theta_{s} * T_{z} ) + E_{down}] * T_{v} / (\pi * d^{2})\]

Таким чином, енергетична світність, відбита від атмосфери, становить:

\[L_{p} = M_{L} * DN_{min} + A_{L} - 0.01* [(ESUN_{\lambda} * cos\theta_{s} * T_{z} ) + E_{down}] * T_{v} / (\pi * d^{2})\]

Існує декілька технік DOS (зокрема DOS1, DOS2, DOS3, DOS4), що ґрунтуються на різних припущеннях щодо \(T_{v}\), \(T_{z}\) та \(E_{down}\) . Найпростішим методом є DOS1, який виходить з наступних припущень (Moran et al., 1992):

  • \(T_{v}\) = 1

  • \(T_{z}\) = 1

  • \(E_{down}\) = 0

Таким чином, енергетична світність, відбита від атмосфери, становить:

\[L_{p} = M_{L} * DN_{min} + A_{L} - 0.01 * ESUN_{\lambda} * cos\theta_{s} / (\pi * d^{2})\]

Результуюча відбивальність земної поверхні визначається за:

\[\rho = [\pi * (L_{\lambda} - L_{p}) * d^{2}]/ (ESUN_{\lambda} * cos\theta_{s})\]

В наступній таблиці наведено значення ESUN [W /(m2 * \(\mu m\))] для сенсорів Landsat.

Значення ESUN для каналів Landsat

Канал

Landsat 1 MSS*

Landsat 2 MSS*

Landsat 3 MSS*

Landsat 4 TM*

Landsat 5 TM*

Landsat 7 ETM+**

1

1983

1983

1970

2

1795

1796

1842

3

1539

1536

1547

4

1823

1829

1839

1028

1031

1044

5

1559

1539

1555

219.8

220

225.7

6

1276

1268

1291

7

880.1

886.6

887.9

83.49

83.44

82.06

8

1369

* за Chander, Markham, & Helder (2009)

** за http://landsathandbook.gsfc.nasa.gov/data_prod/prog_sect11_3.html


For Landsat 8, \(ESUN\) can be calculated as (from http://grass.osgeo.org/grass65/manuals/i.landsat.toar.html ):

\[ESUN = (\pi * d^{2}) * RADIANCE\_MAXIMUM / REFLECTANCE\_MAXIMUM\]

де значення RADIANCE_MAXIMUM та REFLECTANCE_MAXIMUM наведені в метаданих знімка.

Значення ESUN [W /(m2 * \(\mu m\))] для сенсора Sentinel-2 (наведені в метаданих знімка) містяться в наступній таблиці.

Значення ESUN для каналів Sentinel-2

Канал

Sentinel-2

1

1913.57

2

1941.63

3

1822.61

4

1512.79

5

1425.56

6

1288.32

7

1163.19

8

1036.39

8A

955.19

9

813.04

10

367.15

11

245.59

12

85.25

ESUN [W /(m2 * \(\mu m\))] values for ASTER sensor are illustrated in the following table (from Finn et al., 2012).

Значення ESUN для каналів ASTER

Канал

ASTER

1

1848

2

1549

3

1114

4

225.4

5

86.63

6

81.85

7

74.85

8

66.49

9

59.85


Приклад порівняння відбивальності TOA, DOS1 скоригованої відбивальності та Landsat Surface Reflectance High Level Data Products (підсупутникове спостереження) наведено на Рисунку Спектральні сигнатури пікселя забудови.

_images/reflectance_graph.jpg

Спектральні сигнатури пікселя забудови

Порівняння відбивальності TOA, DOS1 скоригованої відбивальності та Landsat Surface Reflectance High Level Data Products


4.7. Перерахунок у температуру

В цьому розділі наведено основну інформацію щодо способу перерахунку у яскравісну температуру на супутнику, реалізованому у SCP та оцінки температури земної поверхні.

4.7.1. Перерахунок у яскравісну температуру на супутнику

For thermal bands, the conversion of DN to At-Satellite Brightness Temperature is given by (from https://www.usgs.gov/core-science-systems/nli/landsat/using-usgs-landsat-level-1-data-product ):

\[T_{B} = K_{2} / ln[(K_{1} / L_{\lambda}) + 1]\]

де:

  • \(K_{1}\) = залежна від каналу стала термального перерахунку (Вт/м кв. * стерадіан * \(\mu m\))

  • \(K_{2}\) = залежна від каналу стала термального перерахунку (Кельвін)

та \(L_{\lambda}\) спектральна густина енергетичної світності на апертурі сенсора, що вимірюється у Вт/(м кв. * стерадіан * \(\mu m\)).

Сталі \(K_{1}\) та \(K_{2}\) для сенсорів Landsat наведено в наступній таблиці.

Сталі термального перерахунку для Landsat

Стала

Landsat 4*

Landsat 5*

Landsat 7**

\(K_{1}\)

671.62

607.76

666.09

\(K_{2}\)

1284.30

1260.56

1282.71

* за Chander & Markham (2003)

** за NASA (2011)


Для Landsat 8 значення \(K_{1}\) та \(K_{2}\) наводяться в файлі метаданих знімка.


\(K_{1}\) та \(K_{2}\) розраховуються як (Jimenez-Munoz & Sobrino, 2010):

\[K_{1} = c_{1} / \lambda^{5}\]
\[K_{2} = c_{2} / \lambda\]

де (Mohr, Newell, & Taylor, 2015):

  • \(c_{1}\) = перша стала випромінювання = \(1.191 * 10^{-16} W m^{2} sr^{-1}\)

  • \(c_{2}\) = друга стала випромінювання = \(1.4388 * 10^{-2} m K\)

Таким чином, для каналів ASTER \(K_{1}\) та \(K_{2}\) наведено в наступній таблиці.

Сталі термального перерахунку для ASTER

Стала

Канал 10

Канал 11

Канал 12

Канал 13

Канал 14

\(K_{1}\)

\(3.024 * 10^{3}\)

\(2.460 * 10^{3}\)

\(1.909 * 10^{3}\)

\(8.900 * 10^{2}\)

\(6.464 * 10^{2}\)

\(K_{2}\)

\(1.733 * 10^{3}\)

\(1.663 * 10^{3}\)

\(1.581 * 10^{3}\)

\(1.357 * 10^{3}\)

\(1.273 * 10^{3}\)

4.7.2. Оцінювання температури земної поверхні

Several studies have described the estimation of Land Surface Temperature. Land Surface Temperature can be calculated from At-Satellite Brightness Temperature \(T_{B}\) as (Weng et al., 2004):

\[T = T_{B} / [ 1 + (\lambda * T_{B} / c_{2}) * ln(e) ]\]

де:

  • \(\lambda\) = довжина хвилі випроміненої енергетичної світності

  • \(c_{2} = h * c / s = 1.4388 * 10^{-2}\) м K

  • \(h\) = стала Планка = \(6.626 * 10^{-34}\) Дж с

  • \(s\) = стала Больцмана = \(1.38 * 10^{-23}\) Дж/K

  • \(c\) = швидкість світла = \(2.998 * 10^{8}\) м/с

Значення \(\lambda\) для термальних каналів супутників Landsat та ASTER можуть бути розраховані з таблиць в Landsat Satellites та Супутник ASTER.

Several studies used NDVI for the estimation of land surface emissivity (Sobrino et al., 2004); other studies used a land cover classification for the definition of the land surface emissivity of each class (Weng et al. 2004). For instance, the emissivity (\(e\)) values of various land cover types are provided in the following table (from Mallick et al., 2012).

Значення випромінності

Тип земної поверхні

Випромінність e

Ґрунт

0.928

Трава

0.982

Асфальт

0.942

Бетон

0.937


4.8. Перелік посилань

  • Ball, G. H. & Hall, D. J., 1965. ISODATA. A novel method of data analysis and pattern classification. Menlo Park: Stanford Research Institute.

  • Chander, G. & Markham, B., 2003. Revised Landsat-5 TM radiometric calibration procedures and postcalibration dynamic ranges Geoscience and Remote Sensing, IEEE Transactions on, 41, 2674 - 2677

  • Chavez, P. S., 1996. Image-Based Atmospheric Corrections - Revisited and Improved Photogrammetric Engineering and Remote Sensing, [Falls Church, Va.] American Society of Photogrammetry, 62, 1025-1036

  • Congalton, R. and Green, K., 2009. Assessing the Accuracy of Remotely Sensed Data: Principles and Practices. Boca Raton, FL: CRC Press

  • Didan, K.; Barreto Munoz, A.; Solano, R. & Huete, A., 2015. MODIS Vegetation Index User’s Guide. Collection 6, NASA

  • ESA, 2020. Sentinel-1 SAR Definitions. Available at https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/definitions

  • ESA, 2020b. Sentinel-1 SAR Definitions. Available at https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/overview

  • ESA, 2019. A machine learning glossary. Available at https://blogs.esa.int/philab/2019/03/29/a-machine-learning-glossary/

  • ESA, 2015. Sentinel-2 User Handbook. Available at https://sentinels.copernicus.eu/documents/247904/685211/Sentinel-2_User_Handbook

  • ESA, 2013. Sentinel-3 User Handbook. Available at https://sentinels.copernicus.eu/documents/247904/685236/Sentinel-3_User_Handbook

  • Finn, M.P., Reed, M.D, and Yamamoto, K.H., 2012. A Straight Forward Guide for Processing Radiance and Reflectance for EO-1 ALI, Landsat 5 TM, Landsat 7 ETM+, and ASTER. Unpublished Report from USGS/Center of Excellence for Geospatial Information Science, 8 p, http://cegis.usgs.gov/soil_moisture/pdf/A%20Straight%20Forward%20guide%20for%20Processing%20Radiance%20and%20Reflectance_V_24Jul12.pdf

  • Fisher, P. F. and Unwin, D. J., eds., 2005. Representing GIS. Chichester, England: John Wiley & Sons

  • JARS, 1993. Remote Sensing Note. Japan Association on Remote Sensing. Режим доступу: http://www.jars1974.net/pdf/rsnote_e.html

  • Jimenez-Munoz, J. C. & Sobrino, J. A., 2010. A Single-Channel Algorithm for Land-Surface Temperature Retrieval From ASTER Data IEEE Geoscience and Remote Sensing Letters, 7, 176-179

  • Johnson, B. A., Tateishi, R. and Hoan, N. T., 2012. Satellite Image Pansharpening Using a Hybrid Approach for Object-Based Image Analysis ISPRS International Journal of Geo-Information, 1, 228. Available at http://www.mdpi.com/2220-9964/1/3/228)

  • Kruse, F. A., et al., 1993. The Spectral Image Processing System (SIPS) - Interactive Visualization and Analysis of Imaging spectrometer. Data Remote Sensing of Environment

  • Mallick, J.; Singh, C. K.; Shashtri, S.; Rahman, A. & Mukherjee, S., 2012. Land surface emissivity retrieval based on moisture index from LANDSAT TM satellite data over heterogeneous surfaces of Delhi city International Journal of Applied Earth Observation and Geoinformation, 19, 348 - 358

  • Mohr, P. J.; Newell, D. B. & Taylor, B. N., 2015. CODATA Recommended Values of the Fundamental Physical Constants: 2014 National Institute of Standards and Technology, Committee on Data for Science and Technology

  • Moran, M.; Jackson, R.; Slater, P. & Teillet, P., 1992. Evaluation of simplified procedures for retrieval of land surface reflectance factors from satellite sensor output Remote Sensing of Environment, 41, 169-184

  • NASA, 2020. What is Synthetic Aperture Radar?. Available at https://earthdata.nasa.gov/learn/what-is-sar

  • NASA, 2013. Landsat 7 Science Data User’s Handbook. Available at https://www.usgs.gov/land-resources/nli/landsat/landsat-7-data-users-handbook

  • NASA, 2011. Landsat 7 Science Data Users Handbook Landsat Project Science Office at NASA’s Goddard Space Flight Center in Greenbelt, 186 http://landsathandbook.gsfc.nasa.gov/pdfs/Landsat7_Handbook.pdf

  • NOAA, 2020. GOES-R Series. Available at https://www.ncdc.noaa.gov/data-access/satellite-data/goes-r-series-satellites

  • Olofsson, P.; Foody, G. M.; Herold, M.; Stehman, S. V.; Woodcock, C. E. & Wulder, M. A., 2014. Good practices for estimating area and assessing accuracy of land change. Remote Sensing of Environment, 148, 42 – 57

  • Ready, P. and Wintz, P., 1973. Information Extraction, SNR Improvement, and Data Compression in Multispectral Imagery. IEEE Transactions on Communications, 21, 1123-1131

  • Richards, J. A. and Jia, X., 2006. Remote Sensing Digital Image Analysis: An Introduction. Berlin, Germany: Springer

  • Sobrino, J.; Jiménez-Muñoz, J. C. & Paolini, L., 2004. Land surface temperature retrieval from LANDSAT TM 5 Remote Sensing of Environment, Elsevier, 90, 434-440

  • USGS, 2015. Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Level 1 Precision Terrain Corrected Registered At-Sensor Radiance Product (AST_L1T). AST_L1T Product User’s Guide. USGS EROS Data Center.

  • Vermote, E. F.; Roger, J. C. & Ray, J. P., 2015. MODIS Surface Reflectance User’s Guide. Collection 6, NASA

  • Weng, Q.; Lu, D. & Schubring, J., 2004. Estimation of land surface temperature–vegetation abundance relationship for urban heat island studies. Remote Sensing of Environment, Elsevier Science Inc., Box 882 New York NY 10159 USA, 89, 467-483

  • Wikipedia, 2017. k-means clustering. Available at https://en.wikipedia.org/wiki/K-means_clustering