Priloga 1: Končno poročilo Raziskovalni projekt V4-1811 v sklopu Ciljnega raziskovalnega programa »Zagotovimo.si hrano za jutri« KONČNO POROČILO Uporaba podatkov satelitskega sistema Sentinel ter nekateri ostalih podatkov daljinskega zaznavanja za kontrolo neposrednih plačil v kmetijstvu Izvajalec: Kmetijski inštitut Slovenije Poročilo pripravila: dr. Uroš Žibrat in mag. Matej Knapič Ljubljana, 30.10.2020 Kazalo vsebine DS 1 Vodenje projekta 3 DS 2 Analiza uporabnosti različnih prostorskih podatkov in drugih podatkov daljinskega zaznavanja za izboljšanje rezultatov podatkov satelita(ov) Sentinel 4 DS 3 Analiza metod za predobdelavo Sentinel posnetkov in testiranje algoritmov za avtomatsko detekcijo sprememb na kmetijskih zemljiščih in kmetijskih parcelah prijavljenih v ukrepe SKP 6 Metode za predobdelavo Sentinel posnetkov 6 Izbor biofizikalnih indeksov 7 Analiza časovne vrste 8 Avtomatsko zaznavanje košnje in oranja 14 DS 4 Določitev testnih območij in preizkušanje izbranih algoritmov za spremljanje ukrepov SKP 17 Izbor območij 17 Izbor ukrepov SKP 20 Testiranje spektralnih indeksov 21 DS 6 Validacija rezultatov avtomatskih procesov identifikacije sprememb ter izbor uporabnih metod 30 DS 7 Komunikacija z naročnikom oziroma uporabniki 32 DS 1 Vodenje projekta Projekt je v večji meri potekal po načrtu dela. Pojavile so se le težave s predvidenim multispektralnim in hiperspektralnim snemanjem testnih površin. Zaradi resne poškodbe brezpilotnega letalnika in kamere na njem, je bil ta sistem celotno leto 2019 na popravilu in ga nismo mogli uporabljati. Hkrati je operater za letalsko snemanje zamenjal letalo, v katerem smo ugotovili resne težave z ozemljitvijo električne napeljave, zaradi česar ni deloval sistem za pozicioniranje letala (IMU – inertial measurement unit) in smo ga bili prisiljeni zamenjati. Ta sistem je ključnega pomena za delovanje hiperspektralnega sistema, tako da tudi tega nismo mogli izvesti. Zato smo pridobili dostop do podatkov satelitskega sistema Planetscope, z razmeroma visoko prostorsko ločljivostjo. Te podatke smo testirali na istem območju in letih, kot podatke Sentinel-2 in WorldView-3/Geoeye 1. Zaradi pandemije virusa SARS-Covid-19 smo izvedli delavnico o uporabi sistema za spremljanje stanja in avtomatsko zaznavanje sprememb na kmetijskih zemljiščih, ter uporabi spletne platforme za vpogled v rezultate, z uporabo spletne platforme Zoom. Delavnice se je udeležilo 25 slušateljev, delovanje sistema so predstavili mag. Matej Knapič, dr. Grega Milčinski in dr. Uroš Žibrat. Datoteke s predstavitve so s strani KIS posredovali slušateljem. DS 2 Analiza uporabnosti različnih prostorskih podatkov in drugih podatkov daljinskega zaznavanja za izboljšanje rezultatov podatkov satelita(ov) Sentinel Spektralne lastnosti poljin so odvisne od številnih dejavnikov, na primer vrste posevka, atmosferskih razmer, topografije in vidnega kota senzorja. Velik del spektralne variabilnosti poljin lahko zajamemo s posnetki v časovni vrsti, zato so slednji tako pomembni za določanje vrst posevkov ter zaznavanje sprememb. Meritve sprememb v času enega leta omogočajo zaznavanje sprememb znotraj ene vrste posevka (na primer košnja travinj) ali spremembe krovnih rastlin (na primer oranje travinj), kot tudi klasifikacijo vrste posevka za nadzor nad izvajanjem SKP. Podatki v časovni vrsti vsebujejo veliko količino informacij, katere lahko analiziramo s številnimi metodami (Preglednica 1). Osnovni statistični podatki v kombinaciji z naprednimi klasifikacijskimi metodami omogočajo določanje vrste posevka, ter zaznavanje sprememb znotraj poljin. Trende merimo v daljših, večletnih obdobjih. Z njimi lahko zaznamo spremembe na večjih območjih, na primer poškodbe in obnova gozdov po žledolomu. Zaradi nizkih prostorskih ločljivosti satelitskih sistemov se pogosto srečamo s t.i. mešanimi spektralnimi podpisi, kjer je v eni rastrski celici zajeto več objektov. V kmetijstvu posamezne rastrske celice zajamejo spektralne podpise senčnih listov, stebelj, cvetov, plodov in ozadja (i.e. zemlje), ki je lahko deloma senčena in različnih vlažnosti. Mešanim spektralnim podpisom se v kmetijstvu ne moremo nikoli popolnoma izogniti, lahko pa zmanjšamo njihov vpliv z uporabo podatkov s čim višjo prostorsko ločljivostjo. Slednja je odvisna od vrste senzorja oziroma satelita, splošno pravilo pa je, da imajo komercialni sateliti boljše karakteristike (Preglednica 2). Intervali snemanja so tudi ključnega pomena za analize posnetkov v časovni vrsti. Predvsem za zaznavanje sprememb na poljinah potrebujemo podatke s čim višjo časovno ločljivostjo. Slednja se prav tako razlikuje med sateliti, trenutno najvišjo imajo sateliti Sentinel. Prostorsko in spektralno najnatančnejše podatke dobimo z uporabo letalskega hiperspektralnega slikanja, ki je možno tudi v časovni vrsti s poljubno ločljivostjo (odvisno od vremenskih razmer in stroškov). Preglednica 1: Dejavniki pridobljeni iz sprektralnih podatkov v časovni vrsti in kaj lahko z njimi merimo. Zadnji stolpec kaže relativno pomembnost posameznega dejavnika med 1 (najmanj pomembno) in 3 (najbolj pomembno) za določanje vrst posevkov in zaznavanje sprememb na kmetijskih zemljiščih v enem letu. Dejavnik Lahko merimo Relativna pomembnost Statistični podatki . Opisna statistika spektralnih 3 podatkov . Tipično za opis podatkov vezanih na sezonske spremembe in fenološki razvoj . Primer: povprečje, najvišje in najnižje vrednosti Metrike . Opisni podatki časovnih odsekov 3 sprememb . Tipično za celoletne podatke . Primer: trajanje, naklon, hitrost sprememb Trend . Večletni trendi 1 . Intervali so lahko neenakomerni Preglednica 2: Primerjava različnih virov daljinskega zaznavanja na večjih območjih. Sentinel-2 WorldView-2 LandSat 8 Pleiades Planetscope Letalsko hiperspektralno Število pasov 12 8 11 5 4 > 100 Ločljivost 10 – 60 m 0,5 – 2,5 m 15 – 100 0,5 – 2 m 3-4 0,47 m m Časovno ločljivost 3 dni 1 dan 16 dni 1 dan 1 variabilno So posnetki Ne* Da Ne Ne Da Da plačljivi? So na voljo Da Da** Da Da Da Da** arhivski posnetki? Relativna 3 1 2 2 1 2 pomembnost * obdelava posnetkov do nivoja L1C ali L2A je plačljiva, odvisno od dobavitelja posnetkov ** arhivski posnetki so na voljo samo na območjih, kjer je v preteklosti nekdo te posnetke naročil DS 3 Analiza metod za predobdelavo Sentinel posnetkov in testiranje algoritmov za avtomatsko detekcijo sprememb na kmetijskih zemljiščih in kmetijskih parcelah prijavljenih v ukrepe SKP Metode za predobdelavo Sentinel posnetkov Multispepktralni podatki satelitov konstelacije Sentinel se hranijo v petih nivojih, glede na stopnjo obdelave podatkov. Stopnja 0 (Level-0) zajema kompresirane surove spektralne podatke z dodatnimi podatki, potrebnimi za nadaljnjo obdelavo (na primer višina in ephemeride senzorja). Na naslednji stopnji, 1A (Level-1A) so na voljo dekompresirani in georeferencirani podatki. Radiometrično korigirani podatki v enotah radiance na vrhu atmosfere (Top-of-atmosphere radiance) so voljo v nivoju 1B (Level-1B). Na tem nivoju je izvedeno tudi popravljanje okvarjenih rastrskih celic z interpolacijo. Ti trije nivoji podatkov niso na voljo za uporabnike satelitskih posnetkov konstelacije Sentinel. Uporabnikom sta na voljo samo nivoja 1B in 2A (Level-1C in Level­2A). Za izdelavo podatkov nivoja 1C se izvede ortorektifikacija posnetkov, torej projekcija posnetkov v kartografsko geometrijo z uporabo digitalnega modela višin, projekciji UTM/WGS84. Podatke se prilagodi na konstantno ločljivost na tleh, glede na nativno ločljivost senzorja v posameznem spektralnem pasu, 10, 20 in 60 m. Podatki so na voljo kot svetlobni odboj (reflektanca) na vrhu atmosfere za vsako rastrsko celico. Na voljo pa so tudi podatki o skupnem ozonu in skupni vodni pari v stolpcu ter povprečnem tlaku nad morsko gladino. Ti podatki so potrebni za izvedbo atmosferskih korekcij in izračun podatkov nivoja 2A. Slednji je izpeljanka nivoja 1C, v katerem so na voljo podatke svetlobnega odboja na dnu atmosfere. V obeh nivojih je možno uporabiti maske oblakov. Maske oblakov so eden ključnih postopkov predprocesiranje satelitskih posnetkov. Vplivajo na to, katere rastrske celice lahko uporabimo za nadaljnje delo ter na njihov spektralni podpis, kar ima neposreden vpliv na razvoj metod in rezultate analiz, na primer klasifikacij posevkov na poljinah in zaznavanje sprememb. Za maskiranje oblakov so najpogosteje uporabljeni štirje algoritmi: Fmask, Sen2Cor, MAJA in s2cloudless (razvoj tega poteka tudi v okviru projekta H2020 Perceptive Sentinel). Algoritem Fmask je bil razvit na podatkih satelitov Landsat in kasneje prilagojen tudi za Sentinel-2. Evropska vesoljska agencija uporablja algoritem Sen2Cor, namensko razvit za podatke satelita Sentinel-2. Ta dva algoritma uporabljata podatke enkratnega zajema posnetkov. Zajem podatkov v časovni vrsti uporablja algoritem MAJA, in je izmed teh treh najzaneslivejša. Vendar zahteva večjo količino podatkov in tehnično težja za uporabo. Ti trije algoritmi v osnovi delujejo podobno, z uporabo različnih numeričnih pragov na podatkih več spektralnih pasov Sentinel-2. Kompleksnejši pristop k zaznavanju oblakov uporablja algoritem s2cloudless. Ta algoritem prav tako deluje na podatkih enkratnega zajema posnetkov, vendar uporablja odločitvena drevesa (algoritem LightGBM) in poda za vsako rastrsko celico verjetnost, da jo pokriva oblak. Ti algoritmi imajo razmeroma podobne zanesljivosti pri določanju oblakov in ločevanju med oblaki in tlemi. Pomembne razlike nastopijo pri določanu visokih oblakov iz rodu cirus, ter zanesljivosti ločevanja med tlemi in oblaki (Preglednica 3). Preglednica 3: Uspešnosti štirih algoritmov za določanje oblakov in ločevanje med oblaki in tlemi. Povzeto po Zupanc, 2017. Fmask Sen2Cor s2cloudless MAJA Oblaki 89,0 % 97,5 % 99,4 % 99,4 % Cirus 88,3 % 87,7 % 83,8 % 64,1 % Tla 92,8 % 94,3 % 97,8 % 99,7 % Pri določanju oblakov in ločevanju med oblaki in tlemi je dosegel algoritem MAJA najboljše rezultate, vendar je bil pri določanju visokih oblakov najslabši. Algoritma Sen2Cor in s2cloudless imata zelo podobne uspešnosti, v povprečju je drugi boljši za prb. 0,5 %. Pomembna razlika med algoritmi Sen2Cor, s2cloudless in MAJA je tudi potreben čas za izračun mask, saj MAJA potrebuje podatke v časovni vrsti in je tako bolj zamuden. Algoritem s2cloudless je prav tako implementiran v orodju eo-learn, zato smo se odločili, da ga uporabimo za predprocesiranje podatkov v tem projektu. Izbor biofizikalnih indeksov Analiz večspektralnih satelitskih posnetkov ne opravljamo na surovih podatkih (svetlobnem odboju, i.e. reflektanci) posameznih valovnih dolžin, temveč delamo s spektralnimi indeksi. Slednji so kombinacije dveh ali več spektralnih pasov in so indikatorji dejavnikov, ki nas zanimajo (na primer fiziološki dejavniki rastlin ali vsebnost dušika v tleh). Delimo jih na vegetacijske, geološke, vodne in vezane na od človeka zgrajene strukture. Eden prvih in še vedno popularnih vegetacijskih indeksov je NDVI, ki je v osnovi enostavna transformacija razmerja med rdečim in bližnje­infrardečim spektralnim pasom (normalized difference vegetation index). Danes imamo na voljo večje število spektralnih indeksov, vsem pa je skupno, da so v osnovi enostavna razmerja med manjšim številom spektralnih pasov, razmerje med dvema ali več indeksi ali pa so prilagoditve obstoječih indeksov. Uporabnost indeksov je odvisna od pasov, ki jih upoštevajo pri izračunu. Indeksi, ki uporabljajo rdeč del spektra, merijo spremembe v vsebnosti klorofila, ki v rdečem in modrem delu spektra močno absorbira svetlobo. Če indeksi uporabljajo pasove v zelenem ali rumenem delu spektra, pa upoštevajo druge pigmente (ksantofile, karotenoide, antocianine). Zanesljivost teh indeksov je odvisna od postavitve in spektralne širine pasov. Na primer, za določanje ksantofilov in karotenoidov je primernejši rumeni pas, ki pa ni na voljo na vseh satelitskih sistemih. Zaradi nizkih prostorskih ločljivosti satelitskih sistemov se pogosto srečamo s t.i. mešanimi spektralnimi podpisi, kjer je v eni rastrski celici zajeto več objektov. V kmetijstvu posamezne rastrske celice zajamejo spektralne podpise senčnih listov, stebelj, cvetov, plodov in ozadja (i.e. zemlje), ki je lahko deloma senčena in različnih vlažnosti. Raziskovalci so zato razvili posebne spektralne indekse, ki kompenzirajo vpliv tal na spektralni podpis (soil adjusted vegetation index – SAVI. Izbor spektralnih indeksov je vezan tudi na želeno prostorsko ločljivost podatkov. Zaradi visoke fragmentiranosti in razmeroma majhnih površin slovenskih poljin, smo se omejili na indekse, ki uporabljajo spektralne pasove z največ 10 m prostorsko ločljivostjo. Izmed indeksov, ki izpolnjujejo ta pogoj, smo jih izbrali enajst (Preglednica 4), ki so vezani na oziroma uporabljajo: (1) razvoj rastlin (bujnost rasti ali indeks listne površine -LAI), (2) razvoj rastlin in so prilagojeni za odboj tal, (3) robni rdeči pas, (4) vsebnost klorofila. Preglednica 4: Izbrani spektralni indeksi. Indeks Namembnost Enačba Referenca OSAVI Prilagojen za odboj tal 1.5*((NIR -RED) / (NIR + RED +0.16)) Rondeaux, Genevieve; Steven, Michael; Baret, Frédéric -Optimization of soil-adjusted vegetation indices 1996 NDVI Bujnost rasti (NIR -RED) / (NIR + RED) Rouse, J.W., Jr.; Haas,R.H.; Schell, J.A.; Deering, D.W. -Monitoring the vernal advancement and retrogradation (green wave effect) of natural vegetation 1973 MCARI Vsebnost klorofila (RE740 -RE705) -(0.2 * (RE740 -GREEN) * (RE740 + RE705)) Daughtry, C. S. T.; Walthall, C. L.; Kim, M. S.; de Colstoun, E. Brown; McMurtrey Iii, J. E. -Estimating Corn Leaf Chlorophyll Concentration from Leaf and Canopy Reflectance 2000 GNDVI Bujnost rasti (NIR -GREEN) / (NIR + GREEN) Gitelson, Anatoly A.; Kaufman, Yoram J.; Merzlyak, Mark N. -Use of a green channel in remote sensing of global vegetation from EOS-MODIS 1996 GOSAVI Prilagojen za odboj tal 1.5*((NIR -GREEN) / (NIR GREEN + 0.16)) + Rondeaux, Genevieve; Steven, Michael; Baret, Frédéric -Optimization of soil-adjusted vegetation indices 1996 PSSRA Vsebnost klorofila NIR / RED Blackburn, G. A. -Spectral indices for estimating photosynthetic pigment concentrations: A test using senescent tree leaves 1998 TCARI Vsebnost klorofila 3*(RE705 -RED) -(0.2 * (RE705 -GREEN) * (RE705 / RED)) Haboudane, Driss; Miller, John R.; Tremblay, Nicolas; Zarco-Tejada, Pablo J.; Dextraze, Louise -Integrated narrow-band vegetation indices for prediction of crop chlorophyll content for application to precision agriculture 2002 TVI LAI 0.5 * (120 * (RE740 -GREEN) ­200 * (RED -GREEN)) Rouse, J.W., Jr.; Haas,R.H.; Schell, J.A.; Deering, D.W. -Monitoring the vernal advancement and retrogradation (green wave effect) of natural vegetation 1973 GLI Vsebnost klorofila ((2 * GREEN) -RED -BLUE) /((2 * GREEN + RED + BLUE)) Gobron, N.; Pinty, B.; Verstraete, M. M.; Widlowski, J. L. -Advanced vegetation indices optimized for up-coming sensors: Design, performance, and applications 2000 MTVI1 LAI 1.2 * (1.2 * (NIR -GREEN) -2.5 * (RED -GREEN)) Haboudane, Driss; Miller, John R.; Pattey, Elizabeth; Zarco-Tejada, Pablo J.; Strachan, Ian B. -Hyperspectral vegetation indices and novel algorithms for predicting green LAI of crop canopies: Modeling and validation in the context of precision agriculture 2004 EVI2 Bujnost rasti 2.4 * (NIR + RED) / (NIR + RED + 1) Miura, T.; Yoshioka, H.; Fujiwara, K.; Yamamoto, H., -Inter-Comparison of ASTER and MODIS Surface Reflectance and Vegetation Index Products for Synergistic Applications to Natural Resource Monitoring 2008 Analiza časovne vrste Pri obdelavi časovne vrste posnetkov lahko govorimo o dveh načinih obdelave in sicer analizi kvalitativnih sprememb spektralnega podpisa v času v rastrski celici ter evidentiranja spremembe v časovni vrsti. Pri ugotavljanju sprememb je metoda sorazmerno jasna, saj se ugotavljajo spremembe z računskimi primerjavami satelitskih podob, medtem ko klasifikacija spektralnih podpisov terja bolj kompleksen pristop. Sateliti Sentinel zajemajo signal na območju Slovenije približno vsake tri dni. Vendar je končna količina razpoložljivih podatkov odvisna od več dejavnikov. Med najpomembnejšimi je oblačnost v času zajema podatkov. Pri prevzemu satelitskih podatkov smo postavili mejo 15% oblačnosti na posameznem posnetku. V enem letu sateliti opravijo približno 100 slikanj, primernih za analize pa je pravilno manj kot polovica. Število primernih posnetkov je tudi odvisno od območja Slovenije. Na primer, okolica Bohinja je pogosto prekrita z oblaki, torej je število primernih posnetkov še manjše. Prav tako lahko pride do okvare podatkov pri prenosu od satelita do bazne postaje, torej so luknje v podatkih naključno razporejene, tako v času (datumih) in prostoru (različna območja Slovenije). Da lahko med seboj primerjamo različna območja in leta, smo manjkajoče podatke interpolirali ločeno za vsako poljino. Tako smo pridobili svetlobni odboj in spektralne indekse za vsak dan v letu. Dodatne pomemben dejavnik sta velikost in oblika poljin, saj je od njiju odvisno število razpoložljivih rastrskih celic. Ker smo želeli pridobiti čim bolj kvalitetne podatke, s čim manj mešanih spektralnih podpisov, smo upoštevali smo tiste rastrske celice, ki so bile v celoti znotraj posamezne poljine (Slika 1). Tako smo sicer zmanjšali število razpoložljivih rastrskih celic, vendar smo hkrati povečali kvaliteto podatkov. Slika 1: Primer dolge in ozke poljine, ker sta v celoti zajeti samo dve rastrski celici (rumeni). Vse ostale celice nosijo mešan spektralni signal. Spektralni indeksi so zasnovani tako, da kar najbolj ocenjujejo biofizikalne parametre rastne odeje ob hkratnem zmanjševanju vpliva zunanjih dejavnikov (npr. vpadni kot sončne svetlobe in vidni kot), normalizirajo interne učinke (variabilnost znotraj rastlinske odeje), ter so vezani na biofizikalne in biološke dejavnike (npr. listna površina – LAI ali vsebnost klorofila). Vendar so vrednosti indeksov med različnimi objekti lahko zelo podobne (Slika 2), kar izrazito omejuje njihovo uporabnost za zanesljivo klasifikacijo celicam). Po drugi strani analiza časovne vrste posameznih indeksov oblikuje nove možnosti razlikovanja, saj so spremembe vrednosti indeksov pogojene z rastjo rastlin (Slika 3). Vendar so metode ločevanja posevkov v časovni vrsti pogoste omejene z uporabo samo enega indeksa, torej so omejene na informacije iz omejenega nabora spektralnih pasov (npr. pri najpogosteje uporabljenem indeksu NDVI na samo dva). Zato smo z analizo glavnih komponent (PCA) združili več indeksov in jih analizirali v časovni vrsti. Klasične multivariantne statistične metode klasifikacij posnetkov obravnavajo posamezne korake ločeno, torej zanemarijo časovno komponento. Za določanje vrste posevkov smo zato testirali več algoritmov strojnega učenja (na primer odločitvena drevesa in nevronske mreže), na podatkih ločenih po posameznih indeksih ter združenih podatkih izbranih indeksov. Slika 2: Primerjava povprečnih vrednosti indeksa MCARI s standardnimi deviacijami v časovni vrsti za vseh šest vrst posevkov za območje Lenart v letu 2017. Slika 3: Povprečja izbranih spektralnih indeksov v časovni vrsti za vseh šest vrst posevkov za območje Lenart v letu 2017. Pri razvoju klasifikacijskih modelov je ključnega pomena zanesljivost oziroma kvaliteta podatkov. Spektralne podatke s satelita je potrebno združiti z bazo podatkov zbirnih vlog. Nadzor nad kvaliteto satelitskih podatkov izvajamo s postopki pred-procesiranja, pri bazi zbirnih vlog pa je nadzor nad kvaliteto omejen. Tako smo pri preliminarnih analizah spektralnih indeksov v časovni vrsti ugotovili, da so v nekaterih zbirnih vlogah navedeni napačni posevki (Slika 4). Najzanesljivejše podatke bi dobili s t.i. ground-truth nadzorom podatkov, kar pomeni, da bi za vse poljine, ki smo jih uporabili za razvoj in validacijo modelov, in vivo preverili stanje (na primer vrsta posevka in aktivnosti na poljini). Na ta način bi imeli zanesljive podatke, s katerimi bi z večjo gotovostjo opravili analize v časovni vrsti. Prav tako bi zmanjšali napake, ki so prisotne v zbirnih vlogah. Zato smo dodali korak predprocesiranja podatkov poljin, kjer so na zbirnih vlogah navedena žita ali koruza. Na teh podatkih smo izvedli nadzorovano klasifikacijo z algoritmom RandomForest. Algoritem je na teh podatkih dosegel 97% pravilnih klasifikacij, pri ločevanju koruze in žit. Razlike v časovni vrsti med koruzo in žiti so razmeroma izrazite, torej je možno natančno ločevanje. Ker smo za vsako poljino dobili verjetnost, da je na njej posamezna vrsta posevka, smo v naslednjem koraku primerjali poljine, kjer sta se razlikovala podatka v zbirnih vlogah in modelna napoved. Če je bila v takšnih primerih verjetnost napačnega posevka visoka (na primer v zbirni vlogi navedena pšenica, model pa je določil koruzo), t.j. nad 90%, smo tej poljini zamenjali vrsto posevka v zbirni vlogi. Tako popravljene podatke smo nato uporabili kot vhodne podatke za razvoj končnih modelov. Trajnih travinj, travnikov na njivi in metuljnic pri tem koraku nismo vključili, saj lahko imajo izrazito pestre vrstne sestave in so spektralno lahko zelo podobne koruzi. Če bi jih vključili v ta korak predpriprave podatkov, bi lahko v podatke vnesli nove napake in tako zmanjšali zanesljivost razvitih modelov. Slika 4: Vrednosti indeksa MCARI v časovni vrsti za območje Lenart v letu 2017 za pšenico. Zgornja slika prikazuje neagregrirane podatke za vse poljine. Spodnja slika prikazuje povprečne vrednosti podatkov združenih po vrsti posevka (pšenica, koruza in neznano). Barve poligonov na spodnji sliki označujejo: modra – pšenica, oranžna – koruza, rumena – neznano. Poleg spektralnih podatkov v obliki spektralni indeksov v časovni vrsti, smo pri razvoju klasifikacijskih modelov uporabili tudi geografske in fenološke podatke. Fenološki razvoj rastlin poteka po dobro opisanih stopnjah (npr. izraščanje listov, hitra rast, stabilna površina), časovna vrsta (hitrost rasti in senescence) tega razvoja pa se razlikuje med različnimi vrstami rastlin. Variabilnost pojava posameznih fenoloških faz je povečana z različnimi datumi setve ter mikroklimatskimi razmerami. Datumi pojava fenoloških faz posameznih posevkov so zato orientacijske narave, predpostavljamo, da so na večini kmetijskih zemljišč določen posevek sejali v časovnem razponi 14 dni. Prav tako je potrebno upoštevati razlike med območji Slovenije (Preglednica 5). Fenološki podatki sami torej niso dovolj zanesljivi za ločevanje posevkov (Slika 5), v kombinaciji s spektralnimi indeksi v časovni vrsti pa lahko izboljšajo klasifikacije vrst posevkov. Na primer, indeks MCARI v primerjavi z NDVI bolj primeren za določanje metličenja koruze in ločevanje med koruzo, travinji in žiti, ter deloma za določanje cvetenja žit (Slika 6). Preglednica 5: Datumi pojava fenoloških faz koruze v letih 2017, 2018 in 2019 za tri območja Slovenije Leto Fenološka faza Osrednjeslovenska Pomurje Dolenjska 2019 2018 2017 vznik 08.05. 07.05. 22.05. 6.listov 14.06. 10.06. 17.06. 9.listov 29.06. 25.06. 01.07. metličenje 13.07. 04.07. 10.07. vznik 03.05. 26.04. 05.05. 6.listov 08.06. 26.05. 03.06. 9.listov 26.06. 09.06. 17.06. metličenje 10.07. 25.06. 02.07. vznik 09.06. 03.05. 02.06. 6.listov 03.07. 08.06. 02.07. 9.listov 17.07. 29.06. 15.07. metličenje 31.07. 13.07. 27.07. Slika 5: Povprečja in standardne deviacije indeksa NDVI v časovni vrsti za šest vrst posevkov za Slika 6: Povprečja in standardne deviacije indeksa MCARI v časovni vrsti za šest vrst posevkov za celotno Slovenijo. Levo leto 2017, desno leto 2018. Za klasifikacije vrst posevka smo uporabili različne vire podatkov, s čimer smo povečali njihovo zanesljivost. Z uporabo več spektralnih indeksov, ki merijo različne fiziološke dejavnike, v časovni vrsti smo v analizah vključili številne dejavnike, ki jih posamezen indeks ignorira. Zaradi klimatskih razmer v Sloveniji in mikroklimatske pestrosti, smo vključili tudi podatke fenološkega razvoja rastlin, ter geografske podatke (t.j. pripadnost posamezne poljine enemu izmed 12 območij). Slovenijo smo razdelili v 12 območij: Dolenjska, Gorenjska, Goriška, Koroška, Notranjska, Obala_Kras, Osrednjeslovenjska, Posavje, Štajerska, Zasavje, Savinjska in Prekmurje. Klasifikacijske modele smo nato razvili na vseh območjih hkrati, da smo zmanjšali verjetnost prevelikega prileganja modelov (angl. overfitting) podatkov in povišali njihovo posploševalno sposobnost (angl. generalizability). Avtomatsko zaznavanje košnje in oranja Satelitski podatki v časovni vrsti omogočajo zaznavanje sprememb na kmetijskih zemljiščih na večjih območjih. Visoka frekvenca slikanja satelitov Sentinel-2 v kombinaciji s sprejemljivo spektralno ločljivostjo omogoča določanje sprememb na kmetijskih poljinah v celotni časovni vrsti. Torej je poleg neposredne primerjave posnetkov določenega območja iz dveh ali treh datumov možno spremljanje sprememb z visoko časovno ločljivostjo. Slednje je ključnega pomena za zaznavanje sprememb na kmetijskih zemljiščih v rastni sezoni, na primer za zaznavanje košnje in oranja. Za zaznavanje slednjih dejavnosti nas zanimajo kratke spremembe vrednosti indeksov. Te spremembe lahko vrednotimo razmeroma natančno s primerjavo absolutnih vrednosti indeksov od različnih datumih in poiščemo datume, ko je prišlo do dovolj velikih sprememb. Vendar je razpon vrednosti posameznega indeksa odvisen od posevka (ali v primeru travinj vrstne strukture), sorte in rastnih razmer (vpliv bioloških in okoljskih dejavnikov). Metoda primerjave absolutnih vrednosti bi torej bila najnatančnejša, če bi imeli na voljo zanesljive, večletne podatke raznih posevkov za vsako poljino. Prav tako je potrebno upoštevati, da imajo različni satelitski senzorji različno postavljene spektralne pasove (različne širine in različni centri pasov, čeprav lahko imajo enaka imena, npr. Rededge ali NIR), kar bo neposredno vplivalo na absolutne vrednosti indeksov. Morebitna izjema je vegetacijski indeks normaliziranih razlik (NDVI), ki je splošni uporabi za različne namene in razmeroma dobro razumljen. Za NDVI obstaja več lestvic vrednosti, na splošno velja, da negativne vrednosti predstavljajo oblake, sneg ali vodo, pozitivne pa tla in rastline. Za slednje so značilne vrednosti nad 0,2, kar je manj nakazuje tla nepokrita z rastlinami. Ne glede na omenjene omejitve uporabe absolutnih vrednosti, lahko z njimi določimo kratke spremembe na kmetijskih poljinah. Namesto primerjave absolutnih vrednosti in določitve mejnih vrednosti oziroma namesto da številčno ovrednotimo razliko v absolutnih vrednostih med dvema datumoma, ovrednotimo hitrost padca vrednosti. To dosežemo z izračunom linearne regresije med dvema datumoma in primerjamo dobljen naklon regresijske linije z eksperimentalno določeno mejno vrednostjo. Tak pristop je neodvisen tako od absolutnih razlik vrednosti indeksov, kot tudi od vrste indeksa (ali meri klorofil, bujnost rasti, LAI itd.). V satelitskih podatkih v časovni vrsti so vedno prisotna nihanja vrednosti, ki so posledica atmosferskih dejavnikov (npr. tankih slojev oblakov) in so torej neodvisna od dejanskega stanja na poljini. Zato je za določanje hitrih sprememb na poljinah potrebno glajenje podatkov v časovni vrsti, da se izognemu lažnim pozitivnim rezultatom. Za izračun linearnih regresij je prav tako ključnega pomena postavitev obeh skrajnih točk (t.j. datumov), saj bi preblizu postavljeni točki privedlo do lažnih pozitivnih rezultatov, preveč narazen postavljeni pa k lažnim negativnim (Slika 7). Slika 7: Primerjava vpliva razdalje med skrajnima točkama za izračun linearne regresije. Levo: razdalja med skrajnima točkama znaša 10 dni. Desno: razdalja med točkama znaša 40 dni. Izmed izbranih spektralnih indeksov se je za določanje vrste posevkov v časovni vrsti obnesel indeks MCARI, zato smo ga uporabili tudi za izračun regresij za določanje hitrih sprememb. Podatke v časovni vrsti smo zgladili z izračunom drsečega povprečja (širina simetričnega okna 10 dni). Ker nismo imeli na voljo podatkov o datumih košenj travinj, smo izračunali povprečne vrednosti indeksa MCARI za vsa travinja, ločeno po testnih območjih. Nato smo na dobljenih podatkih določili točki z največjim padcem vrednosti indeksa. Optimalni časovni razpon 15 dni smo določili s preizkušanjem vrednosti v razponu 10 -30 dni. S časovnim razponom 15 dni odstranimo časovno prekratke padce vrednosti, ki so verjetno posledica atmosferskih dejavnikov. Hkrati je ta razpon dovolj kratek, da zajame padec vrednosti, pred nadaljevanjem rasti rastlin in višanjem vrednosti indeksa. Iz tako dobljenih enačb linearne regresije smo določili mejni naklon, ki nakazuje hitro spremembo vrednosti indeksa. V obdobju od 1.4. do 15.10. na ta način spremljamo pojave košnje. V celotnem letu spremljamo tudi hkratni pojav hitrih sprememb in padca indeksa NDVI pod vrednost 0,2, kar nakazuje, da na poljini ni prisotnih rastlin, kar interpretiramo kot pojav oranja. Pri pojavi oranja štejemo tudi dneve, ko je vrednost indeksa NDVI ostala med 0 in 0,2. Predpostavljamo, da traja vsaj 30 dni, preden se razvije nova rastlinska odeja in se vrednost indeksa dvigne nad 0,2. Kot pojav oranja štejemo primere, kjer je trajanje nizkih vrednosti indeksa NDVI trajalo vsaj 30 dni. V letih 2017 in 2018 na izbranih območjih na travinjih nismo zaznali poljine, ki ne bi bila košena (Preglednica 6). V obeh letih je bila približno polovica poljin košena 3-krat, večina (t.j. več kot 90%) poljin je bila košena dva do štirikrat (Slika 8). V porazdelitvi števila košenj med letoma nismo našli statistično značilnih razlik (Chi2 = 1941,2; p > 0.05, df = 4). Na večini travinj nismo zaznali oranja (Preglednica 7). Med letoma je sicer razmeroma visoka razlika v odstotku travinj v razredih 0 in <30 (število dni, ko je bil zaznan padec vrednosti indeksa MCARI in so bile vrednosti indeksa NDVI pod 0,2), vendar ta razlika ni statistično značilna (Chi2 = 1943,2; p > 0.05, df = 3). V obeh letih smo zaznali oranje na približno 2,5 odstotkih poljin (Preglednica 7). Preglednica 6: Število košenj v letih 2017 in 2018 za izbrana območja. leto št. košenj št. poljin odstotek 1 170 1.02% 2 4205 25.27% 2017 3 8178 49.15% 4 3660 22.00% 5 426 2.56% 1 105 0.58% 2 2475 13.72% 2018 3 10071 55.84% 4 4949 27.44% 5 434 2.41% Preglednica 7: Oranje na travinjih v letih 2017 in 2018 na izbranih območjih. Razred pomeni število dni, ko so bile vrednosti indeksa NDVI pod mejo 0,2. razred 2017 2018 2017 2018 0 10912 15315 65.58% 84.92% <30 5384 2290 32.36% 12.70% 30 343 429 2.06% 2.38% skupaj 16639 18034 Datum Slika 8: Primer poljine, ki je bila košena dvakrat v letu 2017. Oba dogodka košnje sta označena z rdečim poligonom. Modre točke označujejo datume, ko je bila zaznana košnja in rdeča črta predstavlja linearno regresijo skozi ti dve točki. DS 4 Določitev testnih območij in preizkušanje izbranih algoritmov za spremljanje ukrepov SKP Za testiranje algoritmov za spremljanje ukrepov SKP smo skupaj s predstavniki ARSKTRP in MKGP izbrali šest območij, ki odražajo pestrost pridelave in ukrepov SKP, ter posestno strukturo v Republiki Sloveniji. Pri izboru območij je bila pomembna tudi razpoložljivost podatkov. Teh šest območij odraža stanje v celotni državi, kar bo olajšalo aplikacijo algoritmov tudi na drugih območjih. Vseh šest območij je kvadratne oblike s površino 100 km2 (10x10 km). Podatki so na voljo na nivoju GERK-ov in poljin (vrsta kmetijske rastline, kmetijska dejavnost, agrotehnični ukrepi, zaraščenost itd.). Izbor območij Izbrali smo šest območij (od Vzhoda proti Zahodu): Goričko, Lenart, Brežice, Barje, Bohinj, Kras (Slika 9). Območja so kvadratne oblike s stranicami 10x10 km, površina posameznega kvadranta je torej 100 km2. Slika 9: Izbranih šest testnih območij. Od Vzhoda proti Zahodu si sledijo Goričko, Lenart, Brežice, Barje, Bohinj in Kras. Skupno, v obeh letih, je v kvadrantih 75069 poljin oziroma 68266 GERK-ov. Največ poljin je v kvadrantu Goričko, najmanj v Bohinju (Preglednica 8). Poljin je več kot GERK-ov, saj je lahko več kot ena poljina na istem GERK-u. V Lenartu je povprečna površina poljin največja, najmanjša v Bohinju. Preglednica 8: Število in povprečna površina (ar) GERK-ov in Poljin za vseh šest izbranih kvadrantov. 2017 Število Pov. površina Kvadrant Poljina GERK Poljina GERK Barje 8059 7640 57,0 64,9 Bohinj 3897 3870 39,7 39,9 Brežice 6737 5792 76,6 89,4 Goričko 12169 9623 49,7 58,7 Kras 4822 4599 43,3 45,9 Lenart 6695 5778 90,4 102,8 2018 Število Pov. Površina Poljina GERK Poljina GERK 7353 7149 56,6 62,0 3509 3502 43,1 43,3 5204 4888 77,1 85,9 8225 7423 5301 59,5 3309 3275 45,8 47,4 5090 4727 93,0 102,9 Najpogostejša tipa rabe so njive in trajna travišča (identifikatorja 1100 in 1300). Njiva za rejo polžev (1150) je prisotna samo v kvadrantu Barje, oljčnik (1230) na Krasu, matičnjaki (1212) pa so prisotni samo v Brežicah in na Krasu. Med pogostejšimi rabami so tudi vinogradi (1211) in travniški oziroma ekstenzivni sadovnjaki (1222). Vse ostale vrste rabe so zastopane z manj kot 400 poljinami (Preglednica 9). Preglednica 9: Število poljin glede na vrsto rabe zemljišča za vseh šest kvadrantov za leto 2017. Identifikator rabe zemljišč Kvadrant 1100 1131 1150 1170 1180 1190 1211 1212 1221 1222 1230 1240 1300 1320 1610 BARJE 3718221 0 0 0 1 0 6780 042262 5 BOHINJ 42033 0 0 0 0 0 0 31220 032141050 BREZICE 434928 0 41 5 22141 9975 0 119150 7 GORICKO 823525 0 4 1 21860 67377021930432 8 KRAS 89540 0 0 12 49263 5193 4 126849910 LENART 359749 0 9 13 16187 0 69262 0 1562332 2 3 Vsota 21214 197 1 54 31 24 1514 4 295 1007 4 377 17414 210 33 Zaradi izjemno velikih površine v kvadrantu Brežic največjo skupno povprečno površino zavzemajo kmetijska zemljišča v pripravi (1610), sledijo jim intenzivni sadovnjaki (1221), travinja z razpršenimi neupravičenimi značilnostmi (1320), njive (1100) in trajni travniki (1300). V kvadrantu Lenart so v povprečju največje njive, v kv. Brežice največja kmetijska zemljišča v pripravi in intenzivni sadovnjaki, na Krasu pa z gozdnim drevjem porasla kmetijska zemljišča (Preglednica 10). Preglednica 10: Povprečne površine poljin (ar2) glede na vrsto rabe zemljišča za vseh šest kvadrantov za leto 2017. Identifikator rabe zemljišč Kvadrant 1100 1131 1150 1170 1180 1190 1211 1212 1221 1222 1230 1240 1300 1320 1610 BARJE 56.8 34.7 27 0 0 0 1 0 19.7 14.2 0 0 58.2 105.5 46.8 BOHINJ 9.1 17.8 0 0 0 0 0 0 12.3 15 0 0 43.1 92.5 0 BREZICE 85.4 46.7 0 50 19.6 3 20.8 136 178.5 13.8 0 5 60.2 0 379 GORICKO 56.8 22.8 0 14.8 36 2.5 22.4 0 56.8 21.5 0 3.7 38.9 28.5 161.4 KRAS 25 19.3 0 0 97.8 3.5 45.8 68.3 36.2 15 12.8 1 47.1 118.3 31.2 LENART 114.4 24.3 0 8.1 36.2 4.4 31.8 0 113 26.3 0 6.6 73.1 4 12.3 Vsota 347.6 165.6 27 72.8 189.7 13.4 121.8 204.3 416.5 105.7 12.8 16.3 320.6 348.8 630.7 Skupaj s podatki o prijavah smo prejeli tudi podatke o opravljenih nadzorih. Baza podatkov nadzorov nima istih postavk, kot baza podatkov oddanih prijav. Na primer, v bazi opravljenih nadzorov manjka postavka vrsta poljščine oziroma je le-ta navedena samo v povedni, t.j. opisni, obliki. Te podatke bi bilo potrebno izluščiti iz teh povedi in ustrezno urediti, da bi lahko z njimi opravljali analize. Izbrali smo šest vrst v Sloveniji najpogostejših posevkov, ki brez trajnih nasadov pokrivajo približno 98% kmetijskih površin (Preglednica 11). Skupno je bilo v letih 2017 in 2018 s temi šestimi posevki pokritih 54203 poljin. Med njimi je največ trajnega travinja, sledita mu koruza in metuljnice ali ozimna pšenica, odvisno od leta. Med območji vrst posevka niso enakomerno razporejene, na primer v območju Bohinj izrazito prevladuje trajno travinje, podobno kot na Krasu. Preglednica 11: Število poljin z izbranimi vrstami posevka za vsa območja in za obe leti. Vrsta posevka Barje Bohinj Brežice Goričko Kras Lenart Skupaj Koruza 1283 17 1246 1772 26 1219 5563 Metuljnice 341 3 393 696 85 341 1859 Ozimna pšenica 96 0 525 992 24 410 2047 2017 Ozimni ječmen 184 0 465 783 39 390 1861 Trajno travinje 3166 2449 1605 2308 2162 1851 13541 Travnik na njivi 344 30 113 219 123 256 1085 Skupaj 5414 2499 4347 6770 2459 4467 25956 Koruza 1273 20 1220 1871 30 1402 5816 Metuljnice 1645 1 398 566 97 342 3049 Ozimna pšenica 97 0 614 1038 18 618 2385 2018 Ozimni ječmen 208 0 487 777 43 610 2125 Trajno travinje 3221 2423 1613 2318 2180 1894 13649 Travnik na njivi 376 35 135 262 109 306 1223 Skupaj 6820 2479 4467 6832 2477 5172 28247 Skupaj 12234 4978 8814 13602 4936 9639 54203 Izbor ukrepov SKP S predstavniki MKGP in ARSKTRP smo na sestanku, dne 1.4.2019, določili zbrane ukrepe SKP, ki jih bomo spremljali v tem projektu. Poleg izbranih ukrepov smo se tudi dogovorili za nekaj dodatnih ukrepov, katere bomo testirali na manjšem vzorcu, ki je lahko tudi izven testnih območij, odvisno od razpoložljivosti podatkov. Izbrani ukrepi SKP so: . določitev vrste ali skupin vrst kmetijske rastline, . minimalna aktivnost, . košnja, . oranje, . zimsko varstvo tal, . datumsko predpisana prisotnost rastlin. Določanje vrst ali skupine vrst kmetijskih rastlin (angl. crop type) v časovni vrsti je najpomembnejši dejavnik, saj omogoča spremljanje večjega števila ukrepov oziroma zahtev (na primer minimalno aktivnost, določanje datumsko predpisanih prisotnosti določenih kmetijskih rastlin, določanje pokrovnosti). Nekaterih vrst rastlin ni možno ločiti med seboj z uporabo večspektralnih podatkov, zato jih bomo združili v ustrezne skupine (na primer žita). Ločevanje med strnišči in golo zemljo je sicer problematično, saj na strniščih praviloma ni več dovolj organskega materiala in pri spektralni ločljivosti posentkov Sentinel 2 prevlada signal gole zemlje. Določanje mokrišč in šotišč s satelitskimi posnetki Sentinel 2 je praktično nemogoče, predvsem če so vodna telesa plitka. Vodni stolpec je preplitek in senzorji zaznajo podpise bentoške združbe, torej ne vode. Za zaznavanje vode lahko poleg posnetkov Sentinel 2, uporabimo radarske posnetke satelita Sentinel 1. Košnjo bomo spremljali v časovni vrsti, mokrišča in šotišča pa presegajo okvir tega projekta. Testiranje spektralnih indeksov Z vsemi izbranimi indeksi v časovni vrsti smo izvedli klasifikacije vrst poljine z algoritmom odločitvenih dreves RandomForest. Med posameznimi indeksi smo dosegli najvišjo zanesljivost z indeksom MCARI (Preglednica 12). Slednji se je izkazal tudi za najbolj primernega za spremljanje fenološkega razvoja oziroma za zaznavanje razlik v fenološkem razvoju izbranih rastlin (Slika 10). Z indeksom MCARI smo nato testirali različne algoritme strojnega učenja za klasifikacijo vrst rastlin. Najboljše rezultate smo dosegli z indeksom RandomForest in nevronsko mrežo Multilayer Perceptron (Preglednica 13). Za razvoj klasifikacijskih modelov smo izbrali nevronske mreže, saj so odločitvena drevesa nagnjena k prevelikem prileganju podatkom. Spektralne indekse smo združili v komponentah analize glavnih komponent za vsak dan. Ker so nekateri indeksi med seboj visoko korelirani, smo izračunali medsebojne korelacije in pri končni izbiri indeksov za analizo glavnih komponent upoštevali tudi njihovo uspešnost pri klasifikaciji vrst rastlin. Izbrali smo indekse GLI, PSSRA, MCARI, NDVI, TCARI, EVI2, TVI, GNDVI in OSAVI. Pred izračunom glavnih komponent smo vse indekse standardizirali z z-statistiko. Preglednica 12: Rezultati klasifikacij s spektralnimi indeksi v časovni vrsti z algoritmom RandomForest. TP: true positive rate, FP: false positive rate, ROC: receiver operating curve. Indeks Pravilne klasifikacije TP rate FP rate ROC TVI 72% 0,721 0,056 0,948 TCARI 68,5% 0,685 0,063 0,928 PSSRA 71,77% 0,718 0,057 0,943 OSAVI 71,9% 0,719 0,056 0,942 NDVI 70,1% 0,701 0,059 0,941 MTVI1 70,15% 0,701 0,06 0,938 MCARI 73,76% 0,738 0,052 0,947 GOSAVI 71,81% 0,718 0,056 0,939 GNDVI 73,22% 0,732 0,054 0,946 GLI 72,35% 0,723 0,055 0,944 EVI2 70,9% 0,709 0,058 0,937 Preglednica 13: Uspešnost večih algoritmov strojnega učenja pri ločevanju vrst rastlin z uporabi indeksa MCARI v časovni vrsti. TP: true positive rate, FP: false positive rate, ROC: receiver operating curve. Algoritem Pravilne klasifikacije TP rate FP rate ROC SMO 81,26% 0,813 0,094 0,907 SVM 84% 0,84 0,081 0,88 Dl4jMlp 80,1% 0.801 0.1 0.933 RandomForest 89,6% 0,893 0,054 0,978 RF PCA indeksi 67,7% 0,677 0,638 0,676 RF PCA dnevi 74,41% 0,744 0,144 0,926 MPerceptron 88,9% 0,889 0,056 0,964 Klasifikacijske modele smo razvili z uporabo usmerjene umetne nevronske mreže imenovane večplastni perceptron (angl. multilayer perceptron – MLP). Zanj je značilno več plasti vhodnih vozlišč, povezanih kot usmerjen graf med vhodno in izhodno plastjo. Za usposabljanje omrežja pa uporablja povratno širjenje (angl. backpropagation), ki omogoča, da se program strojnega učenja prilagodi glede na preteklo funkcijo. Pri uporabi MLP je potrebno za dosego najboljših rezultatov optimizirati več nastavitev (na primer število skritih slojev). Žal ni zlatega pravila, katere nastavitve bodo dale najboljše rezultate, zato smo postopek treniranja in validacije modelov večkrat ponovili, da smo našli optimalne nastavitve. Podatke smo razdelili v dve skupini: 66% podatkov smo uporabili za treniranje algoritmov in 34% za validacijo. Najprej smo razvili metode za ločevanje vseh šestih vrst posevkov, kjer smo dosegli povprečno uspešnost 81% (Preglednica 14). Najvišjo uspešnost pravilnih klasifikacij smo dosegli pri trajnem travinju in koruzi (obe nad 0.9), najnižjo pri travnikih na njivi (0.89). Ločevanje ozimne pšenice in ozimnega ječmena med seboj je nezanesljivo, prav tako je nezanesljivo ločevanje med vsem tremi vrstami travinj. Na primer, večina travnikov na njivi je bilo uvrščenih med trajna travinja. Razmeroma veliko poljin s koruzo je bilo uvrščeno med metuljnice. Ta rezultat je bil že nakazan pri primerjavi vrednosti indeksov v časovni vrsti, kjer sta distribuciji vrednosti koruze in metuljnic kazali podoben časovni trend in delno prekrivanju v času metličenja. Ječmen in pšenica prav tako rasteta zelo podobno, na travinjih pa je prisotna visoka diverziteta vrst in je ločevanja tako praktično nemogoče. Zato smo vrste posevkov oziroma rastlin združili v tri nove razrede: (1) koruza, (2) travinje (vse tri skupine travinj) in (3) žita (pšenica in ječmen). S tako združenimi podatki smo zvišali uspešnost klasifikacij na 92,74% (Preglednica 15). Uspešnost klasifikacij je bila najnižja pri žitih (0.84), približno 10% procentov jih je bilo napačno uvrščenih med koruzo in približno 6% med žita. Z uporabo podatkov fenofaz rastlin smo zanesljivost povišali na 94,29% (Preglednica 16). Z porabo fenoloških podatkov rastlin smo zvišal uspešnost ločevanja med žiti in travinji ter koruzo. Ker so podatki povezani s podatki KMG-MID in identifikatorji poljin, lahko rezultate prikažemo v obliki kart klasifikacij (Slika 11, Slika 12, Slika 13). Preglednica 14: Rezultati klasifikacije šestih vrst rastlin. Leva stran Preglednica kaže matriko zmede. TP: true positive rate, FP: false positive rate, ROC: receiver operating curve. Klasificirano kot Dejansko koruza metuljnice pšenica ječmen trajno travinje travnik na njivi koruza 9562 644 43 47 147 17 metuljnice 845 2325 377 290 1095 54 pšenica 741 141 2711 919 159 11 ječmen 495 143 758 2770 197 25 trajno travinje 137 158 53 104 26790 39 travnik na njivi 124 240 46 103 1678 215 Preglednica 15: Rezultati klasifikacij vrst rastlin združenih v tri razrede. TP rate FP rate ROC 0.914 0.054 0.981 0.466 0.027 0.92 0.579 0.026 0.946 0.631 0.029 0.951 0.982 0.122 0.973 0.089 0.003 0.799 Klasificirano kot Dejansko Koruza Travinje Žita Koruza 9569 802 175 0.845 0.977 0.907 0.96 0.911 Travinje 881 33067 680 0.962 0.921 0.955 0.932 0.941 Žita 879 520 7630 0.899 0.969 0.845 0.981 0.934 Sensitivity Specificity Pos Pred Value Neg Pred Value Balanced Accuracy Preglednica 16: Rezultati klasifikacij vrst rastlin združenih v tri razrede ob upoštevanju podatkov fenološkega razvoja rastlin. Klasificirano kot Dejansko Koruza Travinje Žita Koruza 9837 525 184 0.962 0.88 0.914 0.963 0.936 Travinje 744 33204 680 0.971 0.975 0.953 0.9372 0.957 Žita 558 405 8066 0.89 0.902 0.849 0.985 0.966 Sensitivity Specificity Pos Pred Value Neg Pred Value Balanced Accuracy je verjetnost prisotnosti za vse vključene vrste rastlin. Temnejše barve predstavljajo višjo verjetnost. Zgoraj: travinje; Sredina: koruza; Spodaj: žita. DS 5 Uporaba visoko ločljivih multispektralnih satelitskih posnetkov in delno tudi letalskega hiperspektralnega slikanja testnih območij ter njihova analiza Ena glavnih hib satelitov konstelacije Sentinel 2 je razmeroma slaba prostorska ločljivost (najmanj 10 m pri pasov v vidnem in bližnje-infrardečem delu spektra). Prav tako so prisotne luknje v časovnem zaporedju posnetkov, ali zaradi vremenskih razmer, ali težav pri prenosu podatkov. Tovrstne luknje lahko zakrpamo z uporabo podatkov drugih satelitov, na primer komercialnih WorldView 2 in GeoEye. Ta satelita imata boljšo prostorsko ločljivost (v vidnem in NIR delu 2,5 m), torej z njima dobimo za isto poljino več podatkov, ob upoštevanju striktnega načela uporabe smo tistih rastrskih celic, ki so v celoti znotraj poljine. Prostorsko in spektralno še natančnejše podatke bi lahko dobili z letalskim hiperspektralnim snemanjem ali z brezpilotnimi letalniki. Podatki konstelacije Sentinel so prosto na voljo in se jih hrani dalj časa v namenskih podatkovnih zbirkah, torej so na voljo tudi arhivski podatki. Podatki WorldView in GeoEye pa so komercialni in so na voljo samo za območja, kjer jih ne kdo naročil, arhiv za celoten svet ne obstaja. Zaradi te omejitve smo za testiranje drugih virov satelitskih podatkov izbrali območje Kras, kjer je bilo na voljo največ posnetkov WorldView-2 in GeoEye. Te posnetke smo obdelali po enakem postopku, kot posnetke Sentinel (ob upoštevanju lastnosti senzorjev WorldView, na primer za izračun radianc), ter podatke preračunali na svetlobni odboj na dnu atmosfere. Ker je bilo posnetkov WorldView + GeoEye premalo za celotno časovno vrsto, smo jih združili s podatki s podatki Sentinel. Podatke WorldView in GeoEye smo torej uporabili za dopolnjevanje podatkov Sentinel. Zaradi razmeroma visokih cen podatkov WorldView, je to tudi realna uporaba teh posnetkov, saj bi bila nabava komercialnih posnetkov v celotni časovni vrsti, torej za celotno leto, predraga. Zaradi visoke prostorske ločljivosti posnetkov drugih satelitov, smo na isti poljini prejeli 20x več podatkov rastrskih celic, kot ob uporabi Sentinel 2. Ker je na satelitu WorldView nameščen tudi pankromatski senzor, je možno s postopkom pankromatskega ostrenja povišati ločljivost multispektralnih posnetkov na 0,5m. V projektu smo načrtovali tudi uporabo letalskega hiperspektralnega slikanja in multispektralnega slikanja z brezpilotnim letalnikom. S tema dvema sistemoma bi prejeli veliko bolj natančne podatke, s katerimi bi dosegli boljše ločevanje vrst posevkov, kot tudi natančnejše določevanje sprememb na poljinah. Zaradi resne poškodbe brezpilotnega letalnika in kamere na njem, je bil ta sistem celotno leto 2019 na popravilu in ga nismo mogli uporabljati. Hkrati je operater za letalsko snemanje zamenjal letalo, v katerem smo ugotovili resne težave z ozemljitvijo električne napeljave, zaradi česar ni deloval sistem za pozicioniranje letala. Ta sistem je ključnega pomena za delovanje hiperspektralnega sistema, tako da tudi tega nismo mogli izvesti. Po predlogu vsebinskih spremljevalk smo vključili druge satelitske podatke, s katerimi smo nadomestili hiperspektralne. Pridobili smo dostop do posnetkov konstelacije Planetscope, z visoko časovno in prostorsko ločljivostjo, vendar s samo štirimi spektralnimi pasovi. Konstelacijo Planetscope sestavlja 120 majhnih satelitov s podobno osnovno, t.j. pred pankromatskih ostrenjem, prostorsko ločljivostjo, kot WorldView-2/Geoeye. Glavna prednost teh posnetkov za ta projekt je razpoložljivost arhivskih posnetkov za vse dni v letu. Komercialni sateliti torej omogočajo zajem izrazito večjega števila rastrskih celic (Preglednica 17, Slika 14). Glavna hiba satelitov konstelacije Planetscope je njihova majhnost, saj vsak posamezen satelit zajema podatke z razmeroma majhna površine, v primerjavi s Sentinel-2 in WorldView-2/Geoeye. Tako je tudi za majhna območja, na primer 10x10 km, potrebno sestavljanje večih rastrskih slik. Na primer, za območje Kras ni bilo na voljo nobenega posnetka, ki bi območje v celoti zajel. Število posameznih posnetkov za območje je bilo med 2 in 9 za en dan. Vsi posnetki niso bili narejeni ob istem času, ampak razporejeni čez posamezen dan. Točen čas zajema satelitskih podatkov vpliva na spektralni odziv, saj je svetlobni odboj, kot ga meri satelit, odvisen od pozicije sonca in senzorja glede na zemeljsko površje ter atmosferskih razmer, ki se dnevno spreminjajo (na primer zaradi močnejše evapotranspiracije poleti ali fosilnih goriv v ogrevalni sezoni). Natančen vpliv te dodatne časovne komponente še ni dovolj dobro raziskan. Pri prevzemu satelitskih posnetkov smo dodali kriterij števila posameznih posnetkov za eno območjo v enem dnevu. Izbirali smo dneve s čim manj posnetki, največ štirimi za posamezen dan. Ob uporabi podatkov satelitov WorldView-2/GeoEye in Planetscope smo malenkost izboljšali uspešnost klasifikacij, iz 90,8% na 91,73% (Preglednica 18 in 19). S podatki konstelacije Planetscope smo dosegli podobno zanesljivost (91,5%), kot ob hkratni uporabi Sentinel-2 in WorldView-2/Geoeye podatkov (Preglednica 20). Pri klasifikacijah so manjše razlike, ki so verjetno posledica časovne neskladnosti zajema posnetkov Planetscope in razlik v postavitvi spektralnih pasov. Pri zaznavanju sprememb na zemljiščih nismo opazili sprememb med viri podatkov (Preglednica 21 in 22). Vse tri konstelacije nimajo postavljenih spektralnih pasov na popolnoma istih valovnih dolžinah, zato lahko pride do razlik v vrednosti svetlobnega odboja. Vendar rezultati kažejo, da so modeli razviti na podatkih Sentinel-2 dovolj robustni in jih lahko uporabimo na drugih satelitskih senzorjih. Preglednica 17: Primerjava števila rastrskih celic na isti poljini med sateliti Sentinel-2, WorldView-2 in Planetscope. Sentinel-2 WorldView-2 Planetscope Najmanjša vrednost 120 3 Prvi kvartil 5 100 15 Mediana 14 280 42 Povprečje 30,43 608,5 91,3 Tretji kvartil 33 660 99 Največja vrednost 2386 42720 7158 Preglednica 18: Rezultati klasifikacije rastlin na Krasu z uporabo podatkov Sentinel. Klasificirano kot Dejansko Koruza Travinje Žita Koruza 41 1 1 0.854 0.999 0.953 0.999 0.927 Travinje 7 4773 22 0.998 0.812 0.994 0.933 0.905 Žita 0 8 83 0.783 0.998 0.912 0.995 0.891 Sensitivity Specificity Pos Pred Value Neg Pred Value Balanced Accuracy Preglednica 19: Rezultati klasifikacije rastlin na Krasu z uporabo podatkov Sentinel in WorldView + GeoEye. Klasificirano kot Dejansko Koruza Travinje Žita Koruza 42 1 1 0.875 0.999 0.954 0.999 0.937 Travinje 6 4773 20 0.998 0.831 0.995 0.934 0.915 Žita 0 8 85 0.802 0.998 0.914 0.996 0.900 Sensitivity Specificity Pos Pred Value Neg Pred Value Balanced Accuracy Preglednica 20: Rezultati klasifikacije rastlin na Krasu z uporabo podatkov Planetscope. Klasificirano kot Dejansko Koruza Travinje Žita Koruza 41 2 1 0.865 0.999 0.953 0.999 0.933 Travinje 6 4775 18 0.998 0.828 0.996 0.935 0.917 Žita 0 9 84 0.80 0.998 0.913 0.996 0.896 Sensitivity Specificity Pos Pred Value Neg Pred Value Balanced Accuracy Preglednica 21: Primerjava števila košenj na Krasu, ob ločeni uporabi Sentinel podatkov in skupaj s podatki WorldView in GeoEye, ter Planetscope. Sentinel Sentinel + WorldView + GeoEye Planetscope št. košenj št. poljin odstotek št. poljin odstotek št. poljin odstotek 0 75 1.52% 75 1.52% 75 1.52% 1 1376 27.88% 1376 27.88% 1376 27.88% 2 2230 45.18% 2230 45.18% 2230 45.18% 3 1254 25.41% 1254 25.41% 1254 25.41% 4 1 0.02% 1 0.02% 1 0.02% 5 75 1.52% 75 1.52% 75 1.52% Preglednica 22: Primerjava rezultatov zaznavanje primerov oranje na travinjih na Krasu, ob ločeni uporabi Sentinel podatkov in skupaj s podatki WorldView in GeoEye, ter Planetscope. Razred pomeni število dni, ko so bile vrednosti indeksa NDVI pod mejo 0,2. Razred 0: 0 dni, 1: 1-29 dni, 2: več kot 30 dni. Sentinel Sentinel + WorldView + GeoEye Planetscope razred vse poljine travinje vse poljine travinje vse poljine travinje 0 4729 95.81% 4729 95.81% 4729 95.81% 1 195 3.95% 195 3.95% 195 3.95% 2 12 0.24% 12 0.24% 12 0.24% DS 6 Validacija rezultatov avtomatskih procesov identifikacije sprememb ter izbor uporabnih metod Modele za klasifikacijo posevkov in detekcijo sprememb na poljinah smo razvili in testiral na podatkih šestih izbranih območij Slovenije, za leti 2017 in 2018. Ker pa so te metode namenjene uporabi na celotnem območju Slovenije, smo izvedli dodatno validacijo modelov s podatki leta 2019 za celotno državo. Za zajem podatkov smo Slovenijo razdelili na 266 kvadrantov (Slika 15). Satelit Sentinel zajame podatke v pasu širokem več sto kilometrov, zato je potrebno te podatke združiti s podatki o poljinah iz zbirnih vlog in narediti presek med njimi. Tako dobimo podatke samo tistih rastrskih celic, ki so za nas zanimive oziroma jih potrebujemo. Za vsakega izmed teh kvadrantov smo prevzeli vse posnetke, ki so ustrezali zahtevam glede pokritosti z oblaki (največ 15%), za vse kvadrante smo tudi izvedli vse zahtevane postopke predprocesiranja, tako da smo delali s spektralnimi podatki na dnu atmosfere. Pri prevzemu podatkov smo naleteli na več težav, na primer okvarjene produkte oziroma posnetke in naključno manjkajoči posnetki. Večja težava je bilo manjkanje naključnih kvadrantov, za katere kljub večim poskusom in popravljanju kode nismo uspeli zajeti podatkov. Izmed 266 kvadrantov nismo uspeli prevzeti podatkov za 45 kvadrantov, večinoma v zahodnem delu Slovenije (Gorenjska in Notranjska). Şkupno smo zajeli podatke za 428.484 poljin. Podatke smo obdelali po metodah, ki smo jih razvili na podatkih za leti 2017 in 2018. Pri validaciji s podatki leta 2019 je algoritem dosegel uspešnost 91,1% (Preglednica 23). Najnižja uspešnost je bila pri žitih (0,87%), najvišja pri travinju (0,95). Z upoštevanjem fenološkega stanja rastlin smo skupno uspešnost povišali za 2,5% na 93,6%. Ponovno smo zvišali uspešnost klasifikacije koruze in žit, prve za 1,0% in druge za 4% (Preglednica 24). Torej lahko z razmeroma visoko zanesljivostjo trdimo, da je določevanje vrste posevka s hkratno uporabo satelitskih, fenoloških in geografskih podatkov razmeroma zanesljivo. Preglednica 23: Rezultati validacije klasifikacij za leto 2019 brez upoštevanja fenofaz. Klasificirano kot Dejansko Koruza Travinje Žita Koruza 54692 16612 16774 0.8810 0.9453 0.7372 0.9785 0.9131 Travinje 10654 243645 10167 0.9309 0.9612 0.9836 0.8475 0.9460 Žita 12493 17761 46920 0.7930 0.9562 0.7423 0.9667 0.8746 Sensitivity Specificity Pos Pred Value Neg Pred Value Balanced Accuracy Preglednica 24: Rezultati validacije klasifikacij za leto 2019 z upoštevanjem fenofaz. Klasificirano kot Dejansko Koruza Travinje Žita Koruza 55104 3576 3756 0,865 0,980 0,883 0,977 0,923 Travinje 1523 304687 3609 0,977 0,957 0,983 0,941 0,967 Žita 7045 3459 46960 0,864 0,972 0,817 0,980 0,918 Sensitivity Specificity Pos Pred Value Neg Pred Value Balanced Accuracy Testirali smo tudi metodo za detekcijo sprememb na celotnem območju Slovenije. Med tem ko v letih 2017 in 2018 nismo zaznali poljin, ki niso bile košene, je bilo takšnih v letu 2019 3737, kar je 0,9% vseh poljin (Preglednica 25). Večino poljin so kosili dva do štirikrat (91,6%). Oranje v zimskem času smo zaznali pri 17.621 poljinah (4,11%), od teh 3.202 na travinju (1,04%). Te poljine niso bile pokrite z rastlinsko odejo najmanj 30 dni (Preglednica 26). Preglednica 25: Število košenj za celotno Slovenijo v letu 2019. št. košenj št. poljin odstotek 0 3737 0.9% 1 11314 2.6% 2 69216 16.2% 3 174742 40.8% 4 148268 34.6% 5 21151 4.9% 6 56 0.0% Preglednica 26: Oranje na poljinah v letu 2019. Rezultati so razdeljeni v razrede, glede na število dni, ko so bile vrednosti indeksa NDVI pod mejo 0,2. Razred 0: 0 dni, 1: 1-29 dni, 2: več kot 30 dni. Število poljin Odstotek razred vse poljine travinje vse poljine travinje 0 381974 294065 89.15% 95.21% 1 28889 11578 6.74% 3.75% 2 17621 3202 4.11% 1.04% DS 7 Komunikacija z naročnikom oziroma uporabniki V začetku novembra, 9.11.2018, smo na KIS pripravili začetni sestanek med AKTRP, MKGP, KIS in Sinergise. Na sestanku smo predstavili program dela s časovnico projekta in se dogovorili o izboru območij in ukrepov SKP, ki jih bomo v projektu testirali (podrobnosti so v DS4). Hkrati z dogovorom poteka sprotna komunikacija med vsebinskima spremljevalkama in izvajalci projekta. Uroš Žibrat in Matej Knapič sta se udeležila konference MARS (Monitoring Agricultural ResourceS) v Dubrovniku, ki jo je organiziral Joint Research Centre (JRC) Evropske komisije skupaj s hrvaško plačilno agencijo. Na konferenci so potekali pogovori o monitoringu ukrepov skupne kmetijske politike z uporabo metod daljinskega zaznavanja, predvsem z uporabo posnetkov satelitov Sentinel. Projekt smo predstavili predstavnikom ZRSVN na sestanku, dne 20.2.2019. Spletna stran projekta je objavljena na naslovu http://www.kis.si/Ciljni_raziskovalni_programi_CRP/Projekt_V4-1811/ V okviru projekta je bila 26. aprila 2019 organizirana predstavitev sistema Sentinel, platforme Sentinel hub in odprtokodne zasnove za analize podatkov daljinskega zaznavanja oziroma potenciala, ki ga nudi orodje EO-learn, ki so ga zasnovali v podjetju Sinergise d.o.o. Predstavitev je podal Grega Milčinski, Sinergise d.o.o. Predstavitve so se udeležili predstavniki ARSKTRP in MKGP, ki so tvorno sodelovali v diskusiji, ki je sledila predstavitvi. Na mednarodnem kmetijskem sejmu AGRA 2019 je bil organiziran posvet z naslovom »Raba satelitskih posnetkov v kmetijstvu – priložnosti in spremembe«, ki sta ga organizirala MKGP in predstavništvo evropske komisije v Slovenije. Na posvetu je projekt predstavil dr. Uroš Žibrat. Osnoven namen in značilnosti projekta so povzeli različni mediji, med ostalim je bil predstavljen v okviru oddaje Ljudje in zemlja. Prispevek je dostopen na internetnem naslovu https://4d.rtvslo.si/arhiv/ljudje-in-zemlja/174635273. O tem dogodku so poročali tudi v drugih medijih: https://ptujinfo.com/novica/politika-gospodarstvo/foto-satelitski-posnetki-pomagajo­ kmetom-ter-izboljsujejo-njihovo. V začetku leta 2020, 14.januarja, smo na KIS organizirali predstavitev rezultatov projekta. Delavnico o uporabi sistema za določanje posevkov in sprememb ter spletne platforme za prikaz rezultat smo izvedli 26.10.2020 z uporabo spletne platforme ZOOM. Delavnice se je udeležilo 28 interesentov z MKGP in ARSKTRP.