Wielotematyczna platforma dla BPN jako przykład wykorzystania danych satelitarnych

Dane pozyskane z pułapu satelitarnego posiadają ogromny potencjał wykorzystania w monitoringu i zarządzaniu obszarami chronionymi. Na ich podstawie możliwe jest obliczanie różnych parametrów środowiskowych, m.in. wskaźników (NDVI, NBR, etc.), które odzwierciedlają kondycję oraz właściwości całego ekosystemu. Niniejsze opracowanie przedstawia wyniki prac mających na celu dostarczenie różnorodnej informacji przestrzennej dla obszaru Biebrzańskiego Parku Narodowego (BPN) w jednym, spójnym serwisie. W szczególności jest to określenie bilansu cieplnego, kondycji roślinności, zasięgu pożaru, zasięgu terenów podmokłych oraz ilości biomasy. Znajomość tych wskaźników może wspomóc osoby zarządzające BPN w podejmowaniu decyzji związanych z realizacją misji ochrony przyrody. Tereny Parku cechują się obecnością licznych rozlewisk, które są ostoją lęgową dla wielu zagrożonych gatunków ptaków (m.in. uszatki błotnej, rybitwy białoskrzydłej, wodniczki) [11]. Jest to obszar zasługujący pod tym względem na wyróżnienie w całej Europie. Warto podkreślić także obecność na tych terenach rzadkich gatunków roślin (m.in. kosaćca bezlistnego i buławnika czerwonego) [1]. Niestety, systematycznie obniżane zdolności retencyjne, czy zmiany klimatyczne sprzyjają powstawaniu środowiska narażonego na pożary i osuszanie gruntów. Efekty prezentowanej pracy wskazują na możliwości wszechstronnego wykorzystania danych teledetekcyjnych i publikowania wyników przetworzeń oraz analiz w ogólnodostępnych serwisach. To z kolei sprzyja popularyzacji udostępniania danych i informacji o środowisku.

Metodyka

Ogólne założenie budowy serwisu opierało się na jak najszerszym wykorzystaniu otwartych danych satelitarnych do analiz środowiska BPN i prezentacji w jednej wspólnej platformie. Punktem wyjścia były problemy i zjawiska występujące w obszarze Parku i jednocześnie możliwe do analizy na danych satelitarnych.  Zdecydowano się na przedstawienie wybranych zjawisk o charakterze przestrzennym w postaci warstw tematycznych:

1. bilans cieplny (temperatura powierzchni),
2. gatunki roślinności (łąki trzęślicowe),
3. kondycja roślinności,
4. zasięg pożaru,
5. zasięg obszarów podmokłych,
6. ilość biomasy (pokrycie terenu).

Każda z warstw tematycznych była opracowana oddzielnie bazując na przyjętej metodyce i niezbędnych danych pozyskanych z lat 2018-2020 z różnych sensorów (optycznych, termalnych i radarowych) i różnych systemów, w tym Landsat i Sentinel (rys. 1). Wykonanie platformy tematycznej wymagało opracowania spójnej bazy danych w technologii Esri i opublikowania w serwisie ArcGIS Online.

Ponieważ opracowanie każdej z warstw tematycznych wiązało się z oddzielnymi pracami badawczymi, w artykule przedstawiono istotne założenia i osiągnięte wyniki dla każdej z nich. Celem głównym niniejszego opracowania jest zaprezentowanie potencjału jaki niosą dane teledetekcyjne oraz możliwości publikacji powstałych z nich produktów na wspólnej  platformie serwisowej dostępnej online.

Realizacja i wyniki prac eksperymentalnych

1. Bilans cieplny

Analiza bilansu cieplnego dla obszaru Biebrzańskiego Parku Narodowego rozpatrywana była w dwóch aspektach. Pierwszym z nich było zbadanie temperatury powierzchni obszaru objętego pożarem w kwietniu 2020 roku. Drugim natomiast była analiza zmiany temperatury w zależności od typu pokrycia terenu. Do przeprowadzenia analizy wykorzystano sceny pochodzące z satelity Sentinel-3 SLSTR pozyskane w zbliżonych terminach: 10.08.2018, 24.08.2019, 06.08.2020.

Głównym zadaniem instrumentu SLSTR znajdującego się na pokładzie Sentinel-3 jest dostarczenie referencyjnego zestawu danych dotyczących temperatury powierzchni lądu i temperatury powierzchni morza dla celów klimatycznych. Produktem wiodącym omawianej analizy był plik danych LST (ang. Land Surface Temperature), zawierający informacje na temat temperatury powierzchni Ziemi. Na rys. 2 przedstawiono porównanie wartości temperatur bezwzględnych dla obszaru BPN na przestrzeni 3 lat, a szczegółowe zestawienie wartości statystycznych temperatur umieszczono w tabeli nr 1.

Zarówno najniższa, jak i najwyższa temperatura wystąpiła w roku 2019. Dla roku 2018 i 2020 otrzymane wartości statystyczne są zbliżone.

W celu przeprowadzenia niezależnej analizy dla obszaru objętego pożarem, obliczono wartości względnej temperatury powierzchni poprzez wyznaczenie odchylenia wartości temperatury każdego piksela od wartości średniej obrazu. Zestawienie otrzymanych wyników przedstawiono w tabeli nr 2.

Pomimo wystąpienia pożaru w kwietniu 2020 roku, w środkowo-wschodniej części Parku nie odnotowano wzrostu wartości temperatur na obejmowanym obszarze. Mogło to być wynikiem szybkiej regeneracji i odbudowy tamtejszej roślinności. Analizę zmian temperatury w zależności od typu pokrycia przeprowadzono w oparciu o dane CORINE Land Cover z 2018 roku. Analizując dane z trzech lat w podobnym terminie, zauważono, że największą wartością temperatury charakteryzują się obszary zabudowane oraz grunty orne, natomiast najniższą wartością – lasy i roślinność krzewiasta w stanie zmian oraz lasy liściaste.

2. Gatunki roślinności

Celem analizy gatunków roślinności było wyznaczenie zasięgu występowania łąk trzęślicowych na podstawie przeprowadzonej klasyfikacji nadzorowanej obrazu satelitarnego Sentinel-2 (10 kanałów spektralnych o rozdzielczości 10 m i 20 m) pozyskanego w dniu 07.08.2020. W celu poprawy wizualizacji zdjęcia satelitarnego wykonano analizę składowych głównych (ang. PCA). Jest to technika transformacji przestrzeni utworzonej z oryginalnych danych obrazowych do nowej przestrzeni obserwacji, tak aby wariancje w pierwszych współrzędnych były największe (nazwane są głównymi składowymi). Analiza ta powoduje, że spektralne różnice są wyraźniejsze i mogą ujawniać subtelne cechy obrazu niedostrzegalne w oryginalnych danych.

Klasyfikację przeprowadzono metodą maksymalnego prawdopodobieństwa zgodnie z zaadaptowaną metodyką za [8]. Do utworzenia pól treningowych wykorzystano udostępnione przez Park warstwy wektorowe z zasięgiem siedlisk, opracowane na podstawie prac terenowych. W następnych etapach wykonano następujące przetworzenia: reklasyfikacja, zamiana rastra na wektor, usunięcie poligonów o powierzchni pojedynczych pikseli, wypełnienie obszarów oraz wygładzenie krawędzi. Ostatecznym rezultatem była warstwa wektorowa z wyznaczonym zasięgiem łąk trzęślicowych (rys. 4).

Dokładność klasyfikacji klasy “łąki trzęślicowe” wyniosła 65%, co jest umiarkowanie zadowalającym wynikiem. Analiza na podstawie zdjęć satelitarnych może być przydatna do zgrubnego wytypowania obszarów łąk trzęślicowych, jednak do dokładniejszego wyznaczenia ich lokalizacji konieczna jest weryfikacja tych obszarów w terenie lub wykonanie przetworzeń w oparciu o hiperspektralne zdjęcia z niskiego pułapu.

3. Kondycja roślinności

Kolejna warstwa tematyczna dotyczyła kondycji roślinności w Biebrzańskim Parku Narodowym, w szczególności w obszarze pożaru. Istnieje szereg wskaźników, służących badaniu kondycji i wigoru roślinności, wskazanych w literaturze jako skuteczne [3][4]. Należą do nich: GEMI (ang. Global Environment Monitoring Index), EVI (ang. Enhanced Vegetation Index), NDVI (ang. Normalized Difference Vegetation Index) oraz NDVI₇₀₅ (ang. Red Edge Normalized Difference Vegetation Index). NDVI to najbardziej popularny wskaźnik opisujący kondycję roślinności. Cechuje go wrażliwość na wpływ atmosfery oraz mała zmienność przestrzenna, ale jednocześnie łatwa interpretacja wartości. Do policzenia wskaźników posłużyły dwie granule Sentinel-2 z 17.10.2019 oraz 21.09.2020, pobrane dla obszaru pożaru. Zdjęcia cechowały się poziomem przetworzenia L2A uwzględniającym korekcję atmosferyczną, a do obliczeń wykorzystano kanały NIR (B8) oraz Red (B4). NDVI przyjmuje wartości z przedziału <-1;1>. Aby przedstawić w bardziej przystępny sposób stan roślinności, ustanowione zostały przedziały wartości (V), którym została przypisana konkretna klasa na podstawie [12, 6]:

V < -0.1                                   wody i chmury

V >= -0.1 ˄ V < 0.1                  gleby odkryte

V >= 0.1 ˄ V < 0.3                   obumarła roślinność

V >= 0.3 ˄ V < 0.5                   sucha i słaba roślinność

V >= 0.5 ˄ V < 0.8                   zdrowa roślinność

V >= 0.8                                  roślinność o najlepszej kondycji

 

Schemat opracowania warstwy tematycznej został przedstawiony na rys. 5.

Na rys. 6 można zaobserwować, że obszar pożaru w 2019 roku, jak wynika z obliczonego wskaźnika NDVI, zdominowany był przez klasy: “sucha i słaba roślinność” (około 65 % terenu) oraz “zdrowa roślinność” (około 35 %). Natomiast we wrześniu 2020 były to klasy: “roślinność o najlepszej kondycji” (około 50 %) oraz “zdrowa roślinność” (około 47 %). Niewielkie fragmenty pokrywała klasa “sucha i słaba roślinność” i były to głównie łąki wykoszone, dlatego te obszary charakteryzowały się słabszą kondycją.

Na obszarze pożaru nie widać spadku kondycji roślinności, co oznacza, że roślinność stosunkowo szybko się zregenerowała. Zróżnicowanie wartości wskaźników nie jest związane bezpośrednio z pożarem – odnosi się do innych kwestii (m. in. do różnego stanu fenologicznego roślinności w poszczególnych latach). Roślinność zdążyła się odnowić po pożarze – wskaźnik roślinności na ogół był wyższy we wrześniu 2020 niż w październiku poprzedniego roku.

4. Zasięg pożaru

Do wyznaczenia obszarów zdegradowanych przez pożary użyto zdjęć satelitarnych systemu Landsat 8. Tereny z widocznymi skutkami klęski żywiołowej ujawniono korzystając z różnic wskaźników Normalized Burn Ratio (NBR), obliczonych dla dat 20.05.2019 oraz 15.05.2020. Wskaźnik NBR jest bardzo czuły na zmiany w stanie szaty roślinnej oraz w wilgotności i cechach gleby związanych z zaistnieniem pożaru [7], [5]. Dzieje się tak ze względu na użyte do obliczenia jego wartości kanały spektralne (w przypadku Landsat 8 są to kanały B05 i B07).

NBR = NIR – SWIR  / NIR + SWIR

Analizę dokładności wyznaczenia zasięgu pożaru przeprowadzono porównując powierzchnie z warstwy referencyjnej z Copernicus Emergency Management Service [2] opracowanej ze zdjęcia Landsat 8 z 2020 roku (tab. 3). Dokładność wykrywania terenów, na których miał miejsce pożar wyniosła blisko 94%.

Oceniono wizualnie obszar opracowania prezentowany w kompozycji RGB 754. Taka kompozycja pozwala na wytypowanie obszarów ze spaloną roślinnością, która jest widoczna w kolorze jasno-różowym lub czerwonym. Zdrowe lub lekko uszkodzone pożarem obszary z roślinnością uwidaczniają się w barwach zieleni i brązu [9]. Analizy dokonano w obu punktach szeregu czasowego obserwacji. Rys. 7 prezentuje przykład zauważonych różnic pomiędzy rokiem 2019 i 2020.

Na rys. 7 widać wyraźnie zmiany na analizowanym obszarze. Łatwo można dostrzec czerwono-brązowe obszary, które świadczą o wypalonych połaciach roślinności, widocznych jako zdrowe na zdjęciach z roku 2019. Warto także zwrócić uwagę na sytuacyjną zgodność z poligonem referencyjnym (kolor czerwony).

Sprawdzono także korelację wyników z wartościami temperatur obiektów na powierzchni ziemi, pozyskanymi z przetworzeń wartości z kanału 10 sceny satelitarnej. Rys. 8 przedstawia porównanie wartości średnich temperatur obliczone na podstawie zdjęć termalnych na obszarach o różnej intensywności pożaru (klasy od najsłabszej (1) do najintensywniejszej (6), na podstawie [13]).

Wszystkie obliczenia przeprowadzono za pomocą skryptów w języku Python posługując się funkcjami dostępnymi w bibliotece arcpy. Pozwoliło to na automatyzację pracy (rys. 9).

5. Zasięg terenów podmokłych

Kolejnym z rozpatrywanych zagadnień był zasięg obszarów podmokłych na terenie Biebrzańskiego Parku Narodowego. Obszary podmokłe można wyodrębnić przetwarzając dane radarowe i optyczne Sentinel [10]. W tym celu wykorzystano zdjęcie radarowe z systemu Sentinel-1 z 30 marca 2020 r. oraz zdjęcie optyczne z systemu Landsat 8 z 28 marca 2020 r.

Wyznaczenie obszarów podmokłych przeprowadzono w dwóch niezależnych procesach. W pierwszym z nich opracowano maskę wody na podstawie zdjęcia radarowego i przyjętej wartości rozproszenia wstecznego (σ₀ w polaryzacji VV w dB). Przetworzenia, które pozwoliły wyodrębnić tereny podmokłe z tego zdjęcia, wykonano w specjalistycznym oprogramowaniu dedykowanym zdjęciom radarowym. Drugi proces obejmował przeprowadzenie klasyfikacji nadzorowanej zdjęcia optycznego z wykorzystaniem pól treningowych. Zdefiniowano siedem klas, w tym trzy dotyczące terenów podmokłych (woda, rzeki, łąki i mokradła). Koncepcja przetwarzania i opracowania wyników została przedstawiona na rys. 10.

Otrzymane wyniki przetworzeń Sentinel-1 i Sentinel-2 porównano z wybranymi klasami obiektów z Bazy Danych Obiektów Topograficznych (BDOT10k), które przedstawiały: wody powierzchniowe, mokradła i szuwary. Wybrane warstwy posłużyły jako dane referencyjne. Wyniki z obu zdjęć przedstawiały zarówno koryto rzeki, jak i tereny podmokłe. Wskazywały one dużo większy obszar niż klasy z BDOT10k. Badania przeprowadzono dla okresu charakteryzującego się dużą wilgotnością gleb, co może być jedną z przyczyn występujących różnic. Ponadto różnice powoduje przyjęcie uśrednionej granicy wód w BDOT10k. Na rys. 11 przedstawiono wyniki przetworzeń zobrazowań Sentinel-1 oraz Sentinel-2.

Po wyznaczeniu obszarów podmokłych na podstawie obu zdjęć przeprowadzono analizę porównawczą zasięgu terenów podmokłych (przykładowy obszar na rysunku 12). Jak można zauważyć, maski z obu zdjęć w większości pokrywają się i dobrze oddają kształt zbiornika wodnego. W przypadku niektórych obszarów wyniki nie są tak jednoznaczne – różnice pojawiały się na terenach z niską roślinnością, które na zdjęciach radarowych wykazały się wilgotnością, co nie zostało wykryte w klasyfikacji na zdjęciach optycznych.

6. Ilość biomasy

Ostatnia warstwa tematyczna opracowana w ramach projektu przedstawia ilość biomasy. Wykorzystano dane teledetekcyjne pozyskane w zakresie optycznym z Landsat 5 z dnia 10.07.2011 r. Na ich podstawie wyznaczono zasięg lasu młodego i starszego. Następnie na podstawie wartości biomasy na jednostkę powierzchni uzyskanej z literatury [15], obliczono sumaryczną jej wartość dla obszaru BPN.

W oprogramowaniu ArcMap, z wykorzystaniem algorytmu maksymalnego prawdopodobieństwa wykonana została klasyfikacja nadzorowana. Fragment obrazu wejściowego i wynikowego zamieszczono na rys. 13. Wykonana ocena dokładności na niezależnych polach testowych wskazuje, że 79% pikseli zostało przypisanych poprawnie do właściwej klasy pokrycia terenu. Wynik klasyfikacji i wskazane klasy pokrycia posłużyły do obliczenia biomasy wskaźnikiem AGB (tabela 4).

AGB (Above Ground Biomass) jest względną miarą służącą do określenia ilości biomasy w przeliczeniu na jednostkę powierzchni [14]. Zgodnie z wartościami zawartymi w tabeli 4, średnia ilość biomasy (obliczona na podstawie dwóch niezależnych klasyfikacji pokrycia terenu) dla lasów znajdujących się w granicach Biebrzańskiego Parku Narodowego to 1 650 370 t/ha.

7. Baza danych i udostępnienie wyników

Produktami otrzymanymi w wyniku prac jest pięć warstw wektorowych oraz sześć rastrowych, zintegrowanych we wspólnej bazie danych. Bazę udostępniono w serwisie ArcGIS Online jako warstwy internetowe: Feature Layer oraz Tile Layer w konfigurowalnej aplikacji „Interaktywna legenda” – jednej z wielu predefiniowanych aplikacji Esri. Pozwala ona w łatwy sposób ustalać, które treści użytkownik chce w danej chwili widzieć na mapie, a włączenie widoczności danej warstwy uruchamia również panel z objaśnieniem symboli. Przy użyciu dostępnych w tej aplikacji narzędzi, dobrano czytelną symbolizację wartości, aby użytkownik mógł dokładnie przyjrzeć się wynikom na tle podkładowej ortofotomapy. Przykładowe wizualizacje zaprezentowano na rys. 14.

Opracowana platforma i wszystkie opisane warstwy tematyczne dostępne są pod adresem:

https://cipw.maps.arcgis.com/apps/instant/interactivelegend/index.html?appid=2b71d4b3c30840d48711cb3429ea89c4

Podsumowanie

Wielość zaprezentowanych przetworzeń i możliwości analizy danych satelitarnych wskazują na ogromny potencjał ich wykorzystania do badania i oceny środowiska wybranego obszaru. Opracowanie wspólnej platformy, integrującej wyniki analiz jako warstw tematycznych, daje kolejne możliwości porównywania i interpretowania danych. Jest też dobrym sposobem na prezentację i publikację wyników, którymi mogą być zainteresowani zarówno decydenci, jak i społeczeństwo.

Podziękowania

Artykuł powstał w ramach zajęć „Teledetekcyjne źródła danych dla SIP” oraz „Zastosowania SIP” i jest efektem zespołowej pracy studentów specjalności SIP (studia drugiego stopnia, kierunek geodezja i kartografia, edycja 2020/21). Autorzy składają podziękowania pracownikom BPN – Panom Grzegorzowi Kwiatkowskiemu i Krzysztofowi Bachowi za udostepnienie danych i uwagi do opracowanych warstw tematycznych oraz dr hab. inż. K. Osińskiej-Skotak, mgr inż. A. Radeckiej i dr inż. J. Pluto-Kossakowskiej za pomoc w opracowaniu merytorycznym przetworzeń danych satelitarnych i inspirację do przygotowania niniejszej publikacji.

Literatura

1. Biebrzański Park Narodowy. Ochrona zagrożonych gatunków roślin. https://www.biebrza.org.pl/ (dostęp 27.04.2021).
2. Copernicus Emergency Management Service https://emergency.copernicus.eu/ (dostęp 13.11.2020).
3. Dąbrowska Zielińska K. et al., Biophysical Properties of Wetlands Vegetation Retrieved from Satellite Images Institute of Geodesy and Cartography, Remote Sensing Centre, 2004.
4. Evangelides C., Nobajas A., Red Edge Normalised Difference Vegetation Index (NDVI705) from Sentinel-2 imagery to assess post fire regeneration, Remote Sensing Applications Society and Environment, 2020.
5. Holden Z.A, Morgan P., Smith A.M.S., Rollins M; Gessler P.E. Evaluation of novel thermally enhanced spectral indices for mapping fire perimeters and comparisons with fire atlas data. International Journal of Remote Sensing, Vol. 26, No. 21, 2005.
6. Jones G., Vaughan R. A., Remote Sensing of Vegetation. Principles, Techniques, and Applications. Oxford University Press, 2010.
7. Miller J., Thode A., Quantifying burn severity in a heterogeneous landscape with a relative version of the delta Normalized Burn Ratio (dNBR), Remote Sensing of Environment, 2006.
8. Modica G., Pollino M., Solano F., Sentinel-2 Imagery for Mapping Cork Oak (Quercus suber L.) Distribution in Calabria (Italy): Capabilities and Quantitative Estimation, 2019.
9. Perez-Cabello F. et al., Mapping erosion-sensitive areas after wildfires using fieldwork, remote sensing, and geographic information systems techniques on a regional scale, Journal of Geophysical Research: Biogeosciences, 2006.
10. Pluto-Kossakowska J., Osińska Skotak K., Łoś H., Weintrit B., The Concept of SAR Satellite Data Use for Flood Risk Monitoring in Poland, 2017.
11. info. Ostoja Biebrzańska. http://ptaki.info/ (dostęp 27.04.2021).
12. Puliti S., Hauglin M., Breidenbach J., Montesano P., Neigh C.S.R., Rahlf J., Solberg S., Klingenberg T.F, Astrup R., Modelling above-ground biomass stock over Norway using national forest inventory data with ArcticDEM and Sentinel-2 data, 2020.
13. The Landscape Toolbox – Tools & Methods for Effective Land Health Monitoring. Normalized Burn Ratio. https://wiki.landscapetoolbox.org/ (dostęp 13.11.2020).
14. Uykun C., Above-ground biomass and carbon estimations and recommendations for forests in Turkey, 2018.
15. Wójtowicz A., Wójtowicz M., Piekarczyk J., Zastosowanie teledetekcji do monitorowania i oceny produktywności plantacji rzepaku, 2005.