În prezent gospodărirea durabilă a pădurilor a devenit o problemă globală, întrucât omul a înțeles importanța pădurii în menținerea echilibrului… [307715]

1. INTRODUCERE

1.1 Necesitatea cercetărilor

În prezent gospodărirea durabilă a pădurilor a devenit o [anonimizat] a înțeles importanța pădurii în menținerea echilibrului ecosferei. [anonimizat].

Deciziile de gospodărire a [anonimizat], care influențează direct sau indirect dezvoltarea ecosistemelor forestiere.

[anonimizat] (Nicolescu, 2009) precum și tipurile de ecosisteme forestiere (Smith ș.a., 1997). [anonimizat]: stabilitatea arborilor (Bibalani, 2008; Păcurar, 2006), alegerea soluțiilor de exploatare (Ciubotaru, 1998; Borz, 2014), productivitatea mijloacelor de recoltare și colectare (Standgard ș.a., 2014), proiectarea și executarea rețelelor de drumuri forestiere (Bereziuc ș.a., 1987), soluțiile de ameliorare a terenurilor degradate (Ciortuz și Păcurar, 2004), modelarea dezvoltării arboretelor (Pretzsch, 2009) etc.

[anonimizat] ([anonimizat].) [anonimizat] a lua în seamă variabilitatea acestor condiții la nivel de detaliu.

[anonimizat], prin care să se facă trecerea la o „silvicultură pe arbore”. [anonimizat] a măsurilor la nivel de arboret nu mai este suficientă. [anonimizat] „[anonimizat].” (Ciubotaru și Păun, 2014). [anonimizat], prelucrare și reprezentare a [anonimizat].

Geomatica stabilește ca obiectiv cunoașterea la nivel de detaliu a [anonimizat], stocarea, procesarea, analizarea și livrarea datelor geografice. Se reunesc astfel sub o denumire unică metode moderne de topografie ([anonimizat], scanare laser etc.), [anonimizat], [anonimizat], geografie cantitativă ș.a. (Gomarasca, 2009).

Aceste tehnologii moderne sunt în continuă perfecționare și permit poziționarea precisă în spațiu a [anonimizat]: agricultură, [anonimizat], geologie, [anonimizat], hidrologie etc.

[anonimizat] o [anonimizat], de regulă greu accesibile. În 1966 se inițiază proiectul CGIS (Sistemul de Informații Geografice al Canadei) prin care se urmărește cartarea suprafețelor destinate agriculturii și silviculturii (Tomlinson, 1990). Instituțiile implicate în demers recunosc de la bun început necesitatea dezvoltării unor soluții tehnologice noi, metodele manuale fiind inadecvate inventarierii unor suprafețe atât de extinse. Prin acest demers, finalizat în 1971, se dezvoltă o serie de soluții inovatoare care stau la baza evoluției viitoare a tehnologiei: scanarea optică a hărților, conversia raster-vector și invers, sau sistemele de administrare a bazelor de date spațiale (Brimicombe, 2009).

În continuare, geomatica oferă suportul tehnologic pentru rezolvarea a numeroase probleme care fac obiectul silviculturii: cartarea fondului forestier (ca poziție geografică dar și ca structură, caracteristici, stare de sănătate etc.), administrarea pădurilor, inventarierea fondului forestier, urmărirea în timp a evoluției pădurii, întocmirea amenajamentelor, proiectarea infrastructurii sau organizarea lucrărilor de exploatare. Se poate vorbi astfel despre „geomatica forestieră”, ramură de cercetare și aplicare destinată silviculturii, care urmărește „investigarea și cunoașterea pădurilor […] prin achiziția, stocarea, modelarea și manipularea datelor spațiale.” (Boș, 2011).

1.2. Motivația cercetărilor

Obiectul tezei de față este tehnologia LiDAR (Light Detection and Ranging), care presupune cartarea suprafețelor prin scanarea cu impulsuri electromagnetice de tip laser. Dincolo de potențialul general al tehnologiei, aceasta este de un real interes în domeniul forestier. Spre deosebire de fotogrametrie sau teledetecție, impulsurile laser pot pătrunde prin stratul de vegetație până la sol. Se fac deci înregistrări la nivelul superior al coronamentului, pe întreaga structură verticală a pădurii dar și la suprafața terenului. Senzorul LiDAR se montează pe o platformă terestră sau aeriană, alături de un Sistem Inerțial de Măsurare (IMU) și un sistem de poziționare globală prin satelit (GNSS). Posibilitatea instalării ansamblului de măsurare pe platforme aeropurtate (avioane sau elicoptere) înseamnă că nu trebuie să se renunțe la avantajul fotogrametriei sau teledetecției, și anume acoperirea eficientă a suprafețelor extinse. Se preconizează că tehnologia va fi indispensabilă, atât în cercetare cât și în scopuri comerciale, în special în regiunile unde pădurea acoperă areale extinse.

Pe plan internațional, tehnologia LiDAR face deja parte din practica uzuală în domeniul forestier, alături de metodele mai vechi ale geomaticii precum fotogrametria sau teledetecția. La noi în țară, tehnologia este încă în faza de început, în care se fac primele cercetări asupra aplicabilității metodei. Majoritatea acestor cercetări se concentrează asupra potențialului tehnologiei pentru rezolvarea problemelor de inventariere sau estimare a anumitor caracteristici dendometrice. Deși această abordare este una importantă și cu utilitate imediată, considerăm că merită atenție și explorarea tehnologiei LiDAR ca metodă de cartare a suprafeței terestre acoperite de vegetație forestieră. În fapt, înțelegerea și urmărirea evoluției arboretelor nu poate fi completă fără analizarea reliefului în care acestea își desfășoară viața.

Prin obiectivele stabilite, teza se amplasează la convergența între mai multe domenii de cercetare: geomatică, scanare laser și silvicultură. Se urmărește investigarea cât mai amănunțită a potențialului tehnologiei LiDAR de a „ilumina” (la propriu și la figurat) relieful pe care se dezvoltă pădurea, ascuns până acum tehnologiilor fotogrammetrice sau de teledetecție.

În concluzie, putem spune că tehnologia LiDAR se dovedește a fi, alături de celelalte metode ale geomaticii forestiere, un important suport de fundamentare a deciziilor în domeniul forestier. Această tehnologie, fără să înlocuiască deplasarea în teren indispensabilă domeniului, oferă date valoroase, de utilitate teoretică și practică domeniului silviculturii.

1.3. Sumar al cercetărilor desfășurate

Tehnologia scanării laser aeriene (ALS) permite folosirea senzorilor LiDAR la cartarea eficientă a suprafețelor extinse. În urma scanării, rezultă un set de date sub forma unui nor de puncte pentru care se cunosc trei elemente: coordonatele X, Y, și Z într-un sistem de referință.

Transformarea acestui nor de puncte inițial într-un model digital al suprafeței terestre, adecvat scopurilor de cercetare sau practicii, presupune două etape de prelucrare:

a. Filtrarea (sau clasificarea) norului de puncte: operațiune prin care se separă acele puncte aflate la nivelul solului de restul setului de date (adică acele puncte care reprezintă obstacole aflate deasupra terenului – clădiri, vegetație etc.); filtrarea se poate realiza în trei moduri:

i. Automat – caz în care se folosesc algoritmi de filtrare.

ii. Manual – operatorul parcurge o reprezentare tridimensională a norului de puncte și realizează filtrarea pe baza analizei vizuale.

iii. Combinat – rezultatul filtrării automate este îmbunătățit prin corectare manuală.

b. Interpolarea: este procesul prin care norul de puncte filtrat este folosit la generarea unei suprafețe continue, prin care se reprezintă variația coordonatei Z a punctelor (deci un model al suprafeței terestre); interpolarea se realizează exclusiv automat, folosind algoritmi specifici.

Realizarea automată a celor două operațiuni nu se poate face fără erori. În cercetarea de față, aceste erori au fost definite, calculate și analizate din punct de vedere al distribuției statistice și spațiale. S-a analizat de asemenea influența factorilor de teren (înclinare, grad de rugozitate etc.) sau vegetație (gradul de acoperire cu vegetație) asupra variației erorilor.

Pentru procesul de filtrare, erorile au fost estimate în două modalități:

a. Prin compararea rezultatului filtrării automate cu cel al filtrării manuale.

b. Prin compararea rezultatului filtrării automate cu un set de date independent, obținut prin măsurători geo-topografice.

În ceea ce privește interpolarea, erorile au fost determinate prin metoda numită cros-validare: cinci la sută din observații au fost extrase aleatoriu din norul de puncte folosit la interpolare, fiind folosite pentru determinarea acurateții suprafeței interpolate.

2. REPREZENTAREA DIGITALĂ A SUPRAFEȚEI TERESTRE

2.1. Aspecte generale privind stocarea datelor geografice

Reprezentarea sub formă de model a mediului înconjurător nu se poate face fără un grad de abstractizare, prin care se reduce complexitatea realității (Harvey, 2015). În cadrul Sistemelor de Informații Geografice, reprezentarea se face apelând la date stocate în format digital. Acest tip de date au două părți (figura 2.1):

Componenta semantică a obiectelor, reprezentată în general sub formă de atribute organizate tabelar.

Componenta spațială, adică o reprezentare geometrică riguroasă a poziției obiectelor, într-un sistem de poziționare definit.

Fig. 2.1. Legătura într-un Sistem de Informații Geografice dintre

componenta semantică (stânga) și componenta spațială (dreapta).

(sursă: Gomarasca, 2009).

De-a lungul timpului s-au dezvoltat diferite formate de stocare a datelor geografice în format digital, în funcție de scopul pentru care se înregistrează datele, caracteristicile obiectelor reprezentate sau gradul de abstractizare urmărit. Aceste formate (sau modele de stocare) se clasifică în două categorii principale: vector și raster.

2.1.1. Modelul vector

Modelul vector servește la stocarea datelor geografice sub formă de hărți digitale, care descriu un teritoriu folosind următoarele tipuri de elemente geometrice (denumite și primitive):

puncte, prin care se reprezintă adimensional obiecte precum stâlpi, arbori etc.;

linii, definite ca seturi de puncte conectate succesiv, prin care se reprezintă obiectele într-o singură dimensiune (de ex. drumuri, căi ferate, rețele de utilități etc.);

poligoane, adică suprafețe închise delimitate de o linie al cărei punct de final coincide cu cel inițial; sunt folosite la reprezentarea obiectelor în două dimensiuni (de ex. clădiri, parcele etc.);

2.1.2. Modelul raster

Modelul raster modelează suprafețe prin subdivizarea acestora într-o rețea de celule, de regulă cu formă pătrată și dimensiuni egale, organizate sub forma unei matrici (Tamaș și Tereșneu, 2010). Din acest motiv modelele raster se mai numesc și modele rețea.

Un model raster permite reprezentarea unui singur atribut al teritoriului analizat – de exemplu altitudinea, acoperirea suprafeței, categoria de folosință, tipul de sol etc. Fiecare celulă a modelului are asociată o valoare singulară pentru atributul modelat, în timp ce componenta spațială este dată de poziția celulei în cadrul rețelei (număr de rând și coloană).

Principala limitare a modelului raster este că, spre deosebire de modelul vector, nu există o corespondență directă între obiecte și reprezentarea acestora. Spre exemplu, dacă în modelul vectorial o clădire este reprezentată printr-un poligon care delimitează amprenta acesteia, în cazul unei model raster clădirea este reprezentată de un grup de celule învecinate care au atributul „clădire” (figura 2.2).

Modelul raster prezintă însă o serie de avantaje:

este o structură de date simplă, cu stocare facilă și accesare eficientă;

organizarea valorilor pe rânduri și coloane înlesnește procesarea într-un mediu informatic;

modelele pot fi suprapuse cu alte surse de date stocate în format matriceal, precum imaginile satelitare;

editarea valorilor individuale se face cu ușurință, valoarea fiecărei celule fiind independentă de celelalte;

Principala caracteristică a unui model raster este rezoluția spațială (numită și rezoluție geometrică), dată de dimensiunea celulelor (sau suprafața din teren corespondentă unei celule). Rezoluția spațială este principalul element care condiționează nivelul de acuratețe cu care se reprezintă obiectele (cu alte cuvinte, gradul de abstractizare al realității): cu cât dimensiunea celulelor scade, vorbim de o creștere a rezoluției și o reprezentare mai fidelă. Rezoluția este exprimată ca o valoare constantă, de regulă în metri, reprezentând dimensiunea pe axele X și Y a celulelor modelului.

În figura 2.2 se exemplifică modurile diferite în care modelele vector și raster reprezintă realitatea.

Fig. 2.2. Realitatea geografică (a) reprezentată simplificat

folosind un model raster (b), respectiv vector (c).

2.2. Modele folosite la reprezentarea suprafeței terestre

Suprafața terestră este suportul pe care se desfășoară fenomenele geografice. Astfel, reprezentarea reliefului este un domeniu de interes în cadrul Sistemelor de Informații Geografice.

În principiu, integrarea informației topografice într-un sistem informatic necesită două elemente: o sursă de date de altitudine și un set de metode de lucru care să permită modelarea, stocarea, analizarea și interpretarea reprezentărilor.

Particularea modelării suprafeței terestre este faptul că nu sunt necesare reprezentări ale obiectelor geografice bine delimitate în spațiu. În schimb, scopul este surprinderea variației continue a altitudinii terenului. Din acest motiv, pentru modelarea suprafeței terestre se pretează modelele de tip raster. O reprezentare în format digital a suprafeței terenului este denumită Model Digital Altimetric (MDA) sau Model Numeric al Terenului (MNAT), în funcție de caracteristicile acesteia.

2.2.1. Modelul Digital Altimetric (MDA)

Burroughs (2015) definește Modelul Digital Altimetric (MDA) ca o reprezentare regulată sub formă de matrice a variației continue a reliefului în spațiu.

MDA este o structură de tip raster, astfel că fiecarei celule a modelului i se atribuie o valoare unică, reprezentând altitudinea pentru centrul celulei. Variația altitudinii în interiorul suprafeței acoperite de o anumită celulă nu este modelată, celula fiind considerată plană. Putem considera modelul ca reprezentativ doar pentru un set discret de puncte dispuse regulat – adică centrele celulelor (El-Sheimy ș.a., 2005). Rezultă deci că Modelul Digital Altimetric nu este o reprezentare continuă a suprafeței, deși este de regulă perceput astfel. Desigur, atunci când rezoluția spațială a modelului este suficient de mare pentru scopul urmărit, distincția nu are efecte asupra utilizării curente.

2.2.2. Modelul Numeric al Terenului (MNAT)

În plus față de Modelul Digital Altimetric, un Model Numeric al Terenului include nu doar date privind variația altitudinii terenului, ci și o reprezentare explicită a formei suprafeței modelate (Wood, 1996).

Un exemplu de Model Numeric al Terenului este modelul de tip vector denumit TIN (Triangular Irregular Network – Rețea Neregulată de Triunghiuri). După cum îi sugerează și numele, modelul descrie suprafața terenului printr-un set de triunghiuri interconectate, de dimensiuni variabile. Pornind de la un set de date discret pentru care se cunoaște altitudinea (vârfurile triunghiurilor) se modelează o suprafață continuă: altitudinea sau derivatele acesteia (panta, expoziția etc.) pot fi determinate pentru orice poziție din interiorul rețelei de triunghiuri. Diferența dintre un model TIN și MDA, sub aspectul vizualizării, este exemplificată în figura 2.3.

În ceea ce privește obținerea unui Model Numeric al Terenului sub formă de structură raster, aceasta presupune îmbunătățirea unui Model Digital Altimetric prin includerea unor date adiționale. De regulă, acestea sunt reprezentate de vectori care delimitează elemente importante de topografie (Goetz ș.a., 2006), precum direcția schimbărilor abrupte de elevație – văi, creste etc. Includerea acestor elemente este necesară pentru a obține o reprezentare a suprafeței corectă din punct de vedere hidrologic, adică una care să permită modelarea scurgerilor (Quinn et al, 1991).

Facem mențiunea că între MDA și MNAT există un anumit grad de confuzie, definirea termenilor fiind diferită de la un domeniu de cercetare la altul (Podobnikar, 2000). În unele situații termenii sunt chiar considerați sinonimi. În lucrarea de față, pentru evitarea confuziilor, se va folosi exclusiv termenul de Model Digital Altimetric (MDA).

(a) (b)

Fig. 2.3. Reprezentare a suprafeței terestre sub formă de Model Digital Altimetric (a) și Rețea Neregulată de Triunghiuri (b).

2.3. Culegerea datelor necesare modelării suprafeței terestre

Pentru reprezentarea terenului este necesar un set de măsurători, de regulă sub formă de puncte, ale altitudinii. În procesul de modelare a suprafeței terestre culegerea datelor inițiale este etapa cea mai costisitoare și îndelungată. În plus, strategia de alegere a punctelor pentru care se măsoară altitudinea, precum și erorile de determinare a poziției punctelor sunt esențiale pentru calitatea produsului final.

Principalele metode de culegere a datelor privind altitudinea (sau date altimetrice) sunt:

Măsurători terestre

Măsurătorile terestre urmăresc folosirea metodelor și aparaturii specifice geodeziei și topografiei, pentru poziționarea în spațiu a punctelor.

În acest caz, procedura de eșantionare a punctelor poate fi adaptată suprafeței: unde terenul este plat punctele măsurate vor fi mai rare, iar unde complexitatea suprafeței crește se vor determina mai multe puncte. Precizia de măsurare specifică geodeziei sau topografiei este ridicată, astfel că modelele obținute sunt caracterizate de o acuratețe înaltă. Costurile de achiziție a datelor sunt însă ridicate și randamentul este scăzut (Lesky et al, 2002), astfel că metoda este potrivită doar la nivel local, pe suprafețe mici de teren.

Digitizarea hărților analogice

O hartă în format analogic poate conține informații despre relief, de regulă reprezentate sub forma curbelor de nivel. Pentru a transpune datele de altimetrie în format digital se parcurg manual de către un operator curbele de nivel. Operațiunea este foarte lentă, erorile au o probabilitate ridicată de apariție și informația este deseori duplicată sau omisă (Doyle, 1978).

De asemenea, este problematică dispunerea punctelor cunoscute – acestea sunt ordonate de-a lungul curbelor de nivel, astfel că nu există informație despre variația topografiei terenului dintre curbe. Efectul este că variațiile locale de altitudine, care apar între curbele de nivel, nu vor fi modelate.

Fotogrammetrie

Fotogrammetria permite înregistrarea eficientă a datelor pentru suprafețe extinse. Prin procedeul numit stereo-restituție se pot obține stereo-modele care surprind topografia terenului. De pe aceste modele se pot extrage date de altimetrie de către specialiști.

Se pot obține însă doar reprezentări ale suprafeței la nivelul superior al obiectelor care o acoperă – așa-numitele Modele Digitale ale Suprafeței (Digital Surface Models – DSM). Dacă terenul este descoperit, DSM-ul este echivalent unui Model Digital Altimetric (MDA). În situația în care terenul este acoperit de obstacole, între DSM și MDA apar diferențe de altitudine, date de înălțimea obiectelor (figura 2.4).

Fig. 2.4. Diferența între un Model Digital al Suprafeței și un Model Digital Altimetric, în situația terenului acoperit.

Interferometrie radar

Interferometria radar este o metodă activă de teledetecție, care folosește o sursă proprie de energie pentru a culege de la distanță date despre suprafețe. Două antene emit unde radar spre aceeiași suprafață, prin combinarea imaginilor rezultate generându-se un model tridimensional. Înregistrarea datelor se poate face din vehicule aeriene, navete spațiale sau platforme satelitare, astfel că randamentul de culegere a datelor este ridicat. Este metoda prin care s-au cules date pentru modele de elevație cu acoperire aproape globală SRTM. Metoda permite măsurători de precizie, care descriu în detaliu variabilitatea topografiei terenului.

LIDAR

Tehnologia LiDAR (Light Detection and Ranging) presupune folosirea unui emițător de unde laser pentru a determina cu precizie și densitate de eșantionare ridicată poziția în spațiu punctelor. Tehnologia va fi tratată pe larg în capitolele 3-4.

2.4. Produse derivate ale Modelului Digital Altimetric

Folosind ca suport Modelul Digital Altimetric, se pot determina o serie de parametri prin care se caracterizează morfologia suprafeței modelate. Acești parametri se numesc derivate ale MDA (El-Sheimy ș.a., 2005). Pentru determinarea lor, se folosesc algoritmi consacrați care realizează operațiuni de vecinătate – valoarea calculată pentru o anumită celulă va depinde de valorile de altitudine ale celulelor învecinate.

În funcție de produsul folosit la generare, derivatele pot fi de două tipuri:

i. Derivate de ordinul întâi: determinate direct pe baza valorilor de altitudine ale MDA.

ii. Derivate de ordin secund: determinate pe baza unui derivat de ordinul întâi.

În continuare vor fi descrise cele mai comune derivate ale MDA, de interes pentru cercetarea de față. Determinarea produselor derivate se face sub forma unei structuri matriceale de tip raster, suprapusă cu Modelul Digital Altimetric și cu aceeași rezoluție spațială.

a. Derivate de ordinul întâi

i. Înclinarea terenului

Înclinarea este dată de unghiul format de planul natural al terenului cu un plan orizontal imaginar (Strahler, 1973). Cu alte cuvinte, înclinarea reprezintă rata de modificare a altitudinii pentru o anumită deplasare orizontală, exprimată în grade sau procente.

Pornind de la un Model Digital Altimetric cu celule rectangulare, determinarea înclinării terenului se face parcurgând următoarele etape:

se identifică vecinii celulei (8 la număr);

se compară altitudinea celulei curente cu altitudinea celulelor învecinate, pentru a se identifica cea mai mare diferență (figura 2.5);

după identificarea celulei învecinate având cea mai mare diferență de nivel, se calculează panta liniei care unește centrele celor două celule (practic linia de cea mai mare pantă);

Valorile pentru înclinare se vor încadra în intervalul 0-100 (dacă este exprimată în procente), respectiv 0-90 (dacă este exprimată în grade).

ii. Expoziția

Expoziția este dată de linia de intersecție dintre planul înclinat al terenului și un plan orizontal imaginar (Strahler, 1973). Pornind de la un Model Digital Altimetric, determinarea expoziției urmează un procedeu similar cu înclinarea: după ce se identifică celula învecinată față de care celula curentă se află la cea mai mare diferență de nivel, expoziția este dată de orientarea liniei care unește centrele celor două celule (exprimată în grade față de direcția Nord). Expoziția va avea astfel valori în intervalul 0-360 grade.

Conceptual, cele două derivate de ordin întâi orientează în spațiu suprafața terenului – în plan vertical (în cazul pantei) și în plan orizontal (în cazul expoziției).

Fig. 2.5. Exemplu de calcul pentru înclinare: valoarea atribuită celulei centrale se va calcula în funcție de elevația celulei din stânga jos, întrucât față de aceasta se identifică diferența de nivel maximă.

b. Derivate de ordin secund

Curbura în plan. Curbura în profil

Curbura se determină pe baza valorilor de înclinare. Astfel, în timp ce înclinarea descrie rata de modificare a altitudinii, curbura descrie rata de modificare a pantei. Curbura în profil se calculează pe direcția liniei de cea mai mare pantă, în timp ce curbura în plan se calculează pe direcția perpendiculară față de linia de cea mai mare pantă.

Semnul valorii numerice a curburii exprimă tipul acesteia: pozitiv pentru concavitate, respectiv negativ pentru convexitate. Valoarea calculată pentru curbură exprimată gradul de concavitate/convexitate în plan sau profil, cu o valoare nulă reprezentând terenul plan.

2.5. Cercetări în domeniul Modelelor Digitale Altimetrice

Cercetarea în domeniul Modelelor Digitale Altimetrice (MDA) are o istorie îndelungată. Termenul își are originea în cercetările realizate de Charles L. Miller în anii 1950, care urmărea folosirea calculatoarelor la proiectarea de autostrăzipe baza datelor de altimetrie înregistrate prin metode fotogrammetrice (Doyle, 1978). Primele utilizări ale MDA au deci loc în domeniul ingineriei civile (Laflamme, 1959; Gladding, 1964; Rosen ș.a., 1970). Ulterior, modelele altimetrice se răspândesc în domenii tot mai diverse, precum: geomorfologie și geologie (Peucker și Douglas, 1975; Pike, 1988; Shary, 2005; Kalbermatten ș.a., 2012), hidrologie (Mark, 1983; Tarboton, 1991; Casas ș.a., 2006; Niță ș.a., 2016), pedologie (unde se estimează că în jur de 80 la sută din aplicațiile realizate în prezent folosesc MDA pentru cartarea automată a solurilor – Bishop și Minasny, 2005) sau climatologie (Dozier, 1980; Lloyd, 2005; Santos-Alamillos și Vazquez, 2013).

În ceea ce privește cercetarea în domeniul silviculturii, Modelele Digitale Altimetrice pot servi ca date-suport la: organizarea lucrărilor de exploatare (Sessions ș.a., 2006; Raymond, 2012; Đuka ș.a., 2015), modelarea distribuției speciilor (Skidmore, 1989; Torma, 2000) sau planificarea rețelelor de drumuri forestiere (Chung și Sessions, 2001; Iordache ș.a., 2012; Aricak, 2015).

Pe lângă folosirea Modelelor Digitale Altimetrice ca date-suport în diferite domenii de activitate sau cercetare, menționăm alte două direcții de cercetare importante:

folosirea MDA pentru estimarea anumitor parametri prin care se descrie suprafața terestră (Collins, 1973; Dietrich și Wilson, 1993; Lacroix ș.a., 2002; Evans, 2012);

evaluarea și îmbunătățirea gradului de acuratețe al Modelelor Digitale Altimetrice (Ackermann, 1978; Li, 1988; Fisher și Tate, 2006; Chaplot ș.a., 2006; Hohle și Hohle, 2009).

3. ECHIPAMENTE LASER PENTRU DETERMINAREA DISTANȚEI

3.1. Principiul de funcționare al laserului

Laserul este un echipament de emisie a unui fascicul coerent de lumină, caracterizat prin: monocromaticitate, directivitate și intensitate. Din punct de vedere fizic, două unde sunt coerente atunci când au aceeași frecvență și lungime de undă, precum și o diferență de fază constantă. Această proprietate permite focalizarea undelor emise de laser asupra unei suprafețe reduse. De asemenea, coerența permite colimarea undelor (reducerea împrăștierii acestora odată cu îndepărtarea de emițător – directivitatea).

Printre numeroasele aplicații ale tehnologiei se numără măsurarea distanțelor cu un grad ridicat de precizie. Această precizie se poate atinge în principal datorită a două caracteristici ale echipamentelor laser:

Se emit pulsuri de înaltă energie, la intervale de timp relativ scăzute.

Lungimea de undă este redusă, deci se poate atinge un înalt grad de colimare folosind deschideri de dimensiuni mici.

În ultimele decenii, pentru tehnologia de cartare în care se folosesc echipamente laser pentru măsurarea distanței s-a consacrat termenul de LiDAR (Bachman, 1979). Spre deosebire de fotogrammetria clasică (care se folosește de o sursă de energie pasivă – lumina solară), tehnologia LiDAR este activă, adică folosește propria sursă de energie. Această energie este emisă sub formă de unde electromagnetice de tip laser. Astfel înregistrarea datelor se poate face teoretic la orice oră (Baltsavias, 1999).

3.2. Principiul de determinare a distanței

În domeniul măsurătorilor terestre, distanța este elementul care servește, alături de unghiuri, la determinarea poziției în spațiu a punctelor. Cele două procedee de determinare a distanței sunt (Kiss și Chițea, 1997):

procedeul direct: instrumentul de măsurare se aplică direct pe distanța de măsurat;

procedeul indirect: distanța se deduce în funcție de alte elemente, fără a fi parcursă;

Măsurarea indirectă a distanței se face folosind instrumente optice sau electronice. Primele dispozitive electronice de măsurare a distanței folosite în topografie au fost lămpile cu mercur sau tungsten. Începând din anii 1970, acestea au început să fie înlocuite cu echipamente de tip laser (Shan și Toth, 2009).

Pentru a permite măsurarea distanțelor, un echipament laser emite fascicule de lumină. Această emisie se poate face în două moduri:

ca pulsuri de înaltă energie emise individual sau

sub forma unui flux continuu de unde laser.

În funcție de modalitatea de lucru, principiul de măsurare a distanței este diferit, după cum urmează:

Pentru pulsuri individuale, determinarea distanței se face prin măsurarea timpului-de-parcurs (adică timpul scurs între emisia pulsului și recepționarea reflexiei acestuia) cu formula:

[3.1]

unde:

R este distanța între platforma de măsurare și suprafața-obiect;

c viteza luminii;

Astfel, distanța R va fi determinată cu formula:

[3.2]

Rezultă că rezoluția pentru măsurarea distanței (ΔR) este direct proporțională cu rezoluția pentru măsurarea timpului Δt (cea mai mică diferență de timp care poate fi înregistrată):

[3.3]

Pentru emisia continuă de pulsuri, măsurarea distanței necesită modularea intensității energiei emise (de exemplu prin emisia de semnale sinusoidale). În acest caz, timpul-de-parcurs al semnalului se determină pe baza diferenței de fază dintre semnalul emis și cel recepționat:

[3.4]

unde:

f este frecvența semnalului;

Φ diferența de fază;

n numărul de lungimi de undă complete incluse în distanța dintre emițător și suprafața-obiect;

Ignorând efectul acestor n lungimi de undă complete, rezoluția pentru măsurarea distanței (ΔR) va fi:

[3.5.]

unde:

ΔΦ este rezoluția diferenței de fază;

f frecvența semnalului;

În acest caz, se observă că rezoluția pentru măsurarea distanței depinde de un nou parametru: frecvența semnalului emis. Astfel, pentru o rezoluție a diferenței de fază dată, se poate obține o mai bună rezoluție de măsurare a distanței (valori mai mici pentru ΔR) prin creșterea frecvenței semnalului.

Din acest motiv, sistemele care emit fluxuri continue de unde laser și măsoară distanța pe baza diferenței de fază sunt de preferat celor care emit pulsuri individuale.

În practică, folosirea acestora este limitată în principal la echipamentele terestre, unde distanțele de măsurat sunt mai mici. Necesitățile sporite de energie pentru generarea unui semnal continuu îngreunează în prezent folosirea tehnologiei pentru măsurători aeriene sau satelitare, unde distanțele de măsurat sunt relativ mari.

O excepție notabilă este reprezentată de sistemul laser prin flux continuu pentru măsurători aeriene ScaLARS (Scanning Laser Altitude and Reflectance Sensor), dezvoltat în scopuri de cercetare de către Institutul de Navigație al Universității din Stuttgart (Hug, 1994).

3.3. Metode de determinare a distanței

Folosind emițătoare laser, înregistrarea distanțelor se poate face în două modalități: prin profile sau prin scanarea suprafețelor.

Măsurători laser prin profile

Metoda presupune folosirea unui echipament laser pentru a măsura distanțele dintre emitățor și o serie de puncte aflate de-a lungul unei linii (Shan și Toth, 2009). Determinând și unghiurile verticale dintre echipamentul de măsurare și aceste puncte, rezultă profile verticale ale terenului (figura 3.1).

Această metodă de măsurare se poate folosi nu doar de la sol, ci și de pe platforme aeropurtate. În acest caz, emitățorul este fix și îndreptat spre sol, iar deplasarea platformei permite măsurarea distanțelor până la o serie de puncte aflate la sol (așezate de-a lungul direcției de zbor). Dacă se cunoaște localizarea emițătorului la momentul emisiei fiecărui puls (prin cuplarea cu un sistem de poziționare), se poate determina poziția punctelor într-un sistem de referință. Astfel de sisteme care realizează profile aeriene se numesc altimetre laser.

Fig. 3.1. Metoda de măsurare prin profile. Din poziția fixă a senzorului A se măsoară distanțele D1 – D7, care împreună cu unghiurile verticale permit calcularea cotei punctelor B1-B7.

Scanarea laser

Folosind echipamente similare celor pentru măsurători prin profile, dar adăugând un mecanism de scanare (de regulă o oglindă sau prismă rotativă – figura 3.2) se pot carta suprafețe, nu doar profile.

Fig. 3.2. Diferite mecanisme de scanare (din stânga sus, în sensul acelor de ceasornic): două tipuri de oglindă oscilantă, scanner Palmer, scanare prin fibre, poligon rotativ.

(sursa imaginii: Wehr și Lohr, 1999).

În cazul scanării este necesară mișcarea continuă a emitățorului pe două axe:

La echipamentele de măsurare de la sol, mișcarea este realizată de un motor electric care permite rotirea emițătorului în plan orizontal și vertical;

În cazul platformelor aeropurtate, mișcarea de-a lungul uneia dintre axe este asigurată de deplasarea platformei. Pentru mișcarea în plan transversal se folosește un sistem de baleiere. Astfel echipamentul măsoară distanțele până la puncte aflate de-a lungul unor profile orientate transversal față de direcția de zbor (figura 3.3);

Fig. 3.3. Prin mișcarea pe două axe a echipamentului de măsurare rezultă

profile transversale față de direcția de zbor

Tiparul de scanare la sol depinde de topografia terenului, viteza și direcția de zbor precum și de mecanismul de scanare. De exemplu, în cazul unei suprafețe orizontale, oglinzile oscilante produc linii sinusoidale, în timp ce scanarea prin fibre produce linii paralele. Un al treilea tip de scanare folosită în mod uzual este cea cu geometrie eliptică. Tiparul de împrăștiere a punctelor asociat acestor trei mecanisme comune de scanare este prezentat în figura 3.4.

Un alt aspect de menționat este faptul că spațierea dintre punctele aflate pe o linie nu este egală, deoarece viteza unghiulară a mecanismului de scanare nu este constantă.

Fig. 3.4. Tipar la sol pentru trei mecanisme de scanare: oglindă oscilantă fără stabilizare (a), oglindă oscilantă cu stabilizare (b), scanner cu fibre (c), scanner Palmer (d).

(sursa imaginii: Gatziolis, 2008)

În urma scanării rezultă o mulțime de puncte poziționate pe trei axe (X, Y, Z) care acoperă întreaga suprafață scanată. Acest set de date este cunoscut sub numele de nor de puncte (figura 3.5).

Fig. 3.5. Nor de puncte 3D rezultat în urma scanării terenului folosind tehnologia LiDAR

3.4. Clasificarea sistemelor de măsurare cu laser

a. După locul colectării datelor

Măsurătorile laser se pot face de la sol (folosind platforme fixe sau mobile), din aer sau din satelit. Astfel se disting următoarele tehnologii de măsurare:

Scanare terestră de pe platforme fixe – sau TLS (Terrestrial Laser Scanning)

Scanare terestră de pe platforme mobile – sau MLS (Mobile Laser Scanning)

Scanare de pe platforme mobile aeropurtate – sau ALS (Airborne Laser Scanning)

Scanare de pe platforme satelitare – sau SLS (Satellite Laser Scanning)

După modul de măsurare

Pentru simplificarea prezentării principiilor tehnologiei laser, s-a omis o caracteristică importantă a undelor laser: capacitatea de penetrare de către acestea a anumitor obstacle (de exemplu vegetația). Astfel, un anumit puls poate produce multiple reflexii în drumul său spre suprafața terestră. Denumirea uzuală pentru reflexia parțială a semnalului de către un obstacol este de „întoarcere”.

Detecția reflexiilor se face pe baza variației intensității pulsului reflectat. Această intensitate depinde de mai mulți factori (Lefsky et al, 2002):

Puterea totală a pulsului emis

Fracția din amprenta pulsului care se suprapune pe obiectul intersectat

Reflectanța obiectului în lungimea de undă a pulsului laser

Fracția din energia reflectată care este orientată spre emițător

Astfel, interacțiunea dintre pulsul laser și o suprafață cu morfologie complexă (de exemplul coronamentul unui arbore) implică o combinație de energii reflectate la diferite nivele și intensități. Rezultatul este o variabilitate continuă a intensității energiei care ajunge înapoi la echipamentul de măsurare.

În funcție de modalitatea prin care senzorul identifică „întoarcerile”, se disting:

Sisteme de măsurare discretă, care identifică unul sau mai multe puncte de maxim ale curbei de intensitate, pe care le stochează ca reflexii individuale. Istoric, primele sisteme de măsurare laser avea capacitatea de a identifica două întoarceri (sau vârfuri ale curbei de intensitate): primul și ultimul. În ultimele decenii, s-au răspândit sisteme care pot identifica mai mult de 10 reflexii.

Sisteme de măsurare continuă, la care echipamentul digitizează și stochează întreaga funcție de intensitate pentru fiecare puls. Determinarea reflexiilor individuale se va face ulterior, la post-procesarea datelor. În plus, funcția de intensitate poate fi analizată pentru a obține informații suplimentare despre textura și morfologia suprafețelor scanate.

Diferențele conceptuale dintre cele două clase de echipamente sunt ilustrate în figura 3.6.

Fig. 3.6. Modul de colectare a datelor pentru sistemele de măsurare discretă, respectiv continuă.

3.5. Tendințe actuale de cercetare

Domeniul cercetării LiDAR este unul activ, urmărindu-se în continuare perfecționarea tehnologiei din punct de vedere al costurilor, datelor care se pot obține sau al preciziei de măsurare.

Dintre direcțiile curente de cercetare amintim:

Senzorii multispectrali

Senzorii cu fotoni unici (în engl. single-photon)

Senzorii de tip Geiger

În cazul senzorilor LiDAR multispectrali, se emit simultan două sau mai multe pulsuri laser de lungimi de undă diferite. Astfel se măsoară răspunsul spectral detaliat al suprafețelor scanate. Aceste date adiționale pot servi, spre exemplu, la clasificarea categoriei de folosință a terenului (Morsy et al, 2017).

Sistemul LiDAR cu fotoni unici (SPL – Single Photon LiDAR), dezvoltat de compania Sigma Space, urmărește să reducă necesarul de energie pentru măsurători, în special cele aeriene. Dacă la sistemele clasice este necesar un minim de 500-1000 fotoni reflectați spre echipamentul de măsurare pentru ca acesta să poată înregistra semnalul reflectat, sistemul SPL ar putea funcționa cu un minim teoretic de câțiva fotoni. Astfel rezultă nori de puncte de densitate ridicată, iar măsurătorile se pot face de la o altitudine mai mare, deci cu un randament sporit.

Sistemul GML (Geiger-Mode LiDAR) este dezvoltat de compania Harris în scopuri militare. Tehnologia folosește, în locul unui receptor unic, o matrice de componente receptoare de tip diodă (asemănătoare, ca principiu de funcționare, unei camere fotografice digitale). În loc să scaneze gradual suprafața, sistemul o iluminează simultan, capturând informații la o rezoluție ridicat (datorită densității înalte a matricei de diode). Astfel, prin suprapunerea înregistrărilor, obiectele sunt scanate nu doar odată, ci de numeroase ori. Similar sistemului SPL, avantajele teoretice ale tehnologiei sunt o densitate ridicată a punctelor și posibilitatea colectării de date de la o altitudine mai mare, datorită necesarului scăzut de energie.

Se face mențiunea că astfel de sisteme sunt încă în faza de dezvoltare, cercetare și testare, gradul de impact al acestora asupra domeniului practic fiind încă necunoscut.

TEHNOLOGIA DE SCANARE LASER AERIANĂ (ALS)

Tehnologia de scanare laser aeriană (ALS – Airborne Laser Scanning) presupune montarea unui echipament de măsurare a distanței de tip laser pe o platformă aeropurtată (elicopter, avion sau vehicul autonom –„dronă”) pentru cartarea eficientă a suprafețelor extinse.

4.1. Componentele unui sistem ALS

Un sistem tipic de scanare laser aeriană are următoarea structură (figura 4.1) :

Unitate laser de măsurare a distanței, cu sisteme optice pentru transmisia și recepționarea undelor.

Mecanism de scanare (oglindă sau prismă rotativă) care direcționează pulsurile laser de-a lungul unei traiectorii transversale față de direcția de zbor. De asemenea este atașat un dispozitiv care măsoară cu precizie unghiul sub care se emite un puls, denumit encoder unghiular.

Sistem de poziționare globală prin satelit (GNSS – Global Navigational Satellite System) care înregistrează poziția platformei într-un sistem de referință. Aceste date sunt necesare pentru poziționarea punctelor măsurate, dar și pentru a determina elemente necesare la post-procesare, precum altitudinea sau viteza de zbor.

Sistem inerțial de orientare (IMU – Inertial Measurement Unit), care servește la determinarea orientării pe trei axe a platformei de zbor (rotație, ruliu și tangaj).

Unitate hardware/software care permite controlul sistemului, asigură funcționarea corectă a echipamentelor și este responsabilă pentru stocarea datelor pe suport electronic.

Poziția geografică a punctelor scanate (latitudine, longitudine, altitudine) rezultă pe baza următoarelor elemente: i. poziția platformei la momentul emisiei pulsului; ii. înclinarea platformei pe cele trei axe (X, Y, Z) la momentul emisiei pulsului; iii. distanța dintre platformă și punctul care a produs reflexia pulsului laser.

În unele situații, la componentele descrise mai sus se poate atașa un dispozitiv de capturare a imaginilor (cameră fotografică sau video) care oferă date auxiliare.

Fig. 4.1. Componența sistemului de scanare LMS-Q560 dezvoltat de compania Riegl, cu o structură tipică pentru un astfel de echipament.

(după Riegl Technical Brochure, www.riegl.com)

4.2. Surse de erori pentru măsurarea de pe platforme aeropurtate

Sistemele pentru scanarea laser aeriană (ALS) au o structură complexă, cu multiple unități interconectate. Astfel apar mai multe surse de erori, sistematice sau întâmplătoare. Principalele elemente care contribuie la eroarea finală de poziționare a punctelor sunt analizate în continuare:

Senzorul laser

Echipamentului de măsurare i se asociază trei surse de erori principale: eroarea de determinare a distanței, eroarea de determinare a unghiului de scanare și efectul reflectanței suprafeței scanate.

Eroarea de determinare a distanței

Eroarea de determinare a distanței este cea mai dificil de modelat, dar are contribuția cea mai redusă, atunci când sistemul este bine calibrat. Eroarea de măsurare a distanței are efect în principal asupra acurateții de poziționare verticală a punctelor (figura 4.2), iar dacă este sistematică produce deformarea suprafeței scanate (Csanyi, 2007). Detalii tehnice suplimentare despre această sursă de erori se regăsesc în Baltsavias (1999).

Eroarea de determinare a unghiului de scanare

Eroarea de determinare a unghiului de scanare (diferența dintre valoarea măsurată și cea reală) variază sinusoidal cu poziția. Abaterea unghiulară are efect asupra poziționării punctelor în plan vertical, precum și de-a lungul direcției de scanare (figura 4.3). Eroarea poate fi modelată și compensată la calibrarea sistemului LiDAR.

Fig. 4.2. Eroarea de determinare a distanței Δr provoacă o deplasare a

punctului pe axele Y și Z.

Fig. 4.3 Eroarea de determinare a unghiului de scanare Δβ are efect asupra acurateții de determinare a poziției punctului scanat, în plan vertical precum și pe axa Y.

Efectul reflectanței suprafeței scanate

Suprafețele scanate au un răspuns spectral diferit la lungimea de undă a pulsurilor emite de echipamentul de măsurare. Atunci când reflectanța este scazută, reflexia pulsului laser va fi întârziată (distanța estimată va fi mai mare decât cea reală). Apare astfel o eroare verticală negativă: altitudinea calculată a punctului este mai mică decât cea reală. În cazul invers, atunci când suprafața scanată are o reflectanță ridicată, punctele vor fi afectate de o eroare verticală pozitivă.

Întrucât punctele înregistrate au asociată o valoare pentru intensitate (care aproximează, într-o anumită măsură, reflectanța suprafeței), efectul acestui fenomen poate fi diminuat în etapa de procesare a datelor.

Sistemul de poziționare GNSS/IMU

Acuratețea de poziționare a punctelor depinde de două componente ale sistemului de măsurare: echipamentul GNSS (folosit pentru determinarea coordonatelor platformei în sistem global) și echipamentul IMU (folosit pentru determinarea rotației platformei pe cele 3 axe).

Eroarea de poziționare a platformei are efect în principal asupra poziției orizontale a punctelor măsurate. Atunci când scopul colectării datelor este generarea de modele ale suprafeței terestre, capătă importanță creșterea acurateții de poziționare verticală a punctelor, în dauna acurateții de poziționare în plan. Din acest motiv, în trecut accentul în cercetare s-a pus pe modelarea și compensarea celorlalte surse de erori, care afectează poziționarea pe verticală a punctelor. Studiile mai recente atrag însă atenția asupra importanței erorilor de poziționare pe orizontală (Mass, 2001). Dacă efectul acestora este neglijabil în cazul terenurilor plate, situația se schimbă în cazul terenurilor înclinate. În acest caz, eroarea de poziționare orizontală va influența precizia de determinare a altitudinii.

Erori asociate componentei GNSS

Echipamentul GNSS (Global Navigational Satellite System) permite poziționarea platformei de măsurare, și deci a punctelor înregistrate, într-un sistem de poziționare global. Pentru a determina poziția, echipamentul se află în comunicare continuă cu o serie de sateliți.

Dintre factorii care influențează așa-numita „soluție de navigare” oferită de echipamentul GNSS, amintim:

Caracteristicile constelației satelitare (numărul de sateliți și dispunerea acestora pe bolta cerească).

Eroarea satelitară (mici deviații ale sateliților de la orbitele calculate).

Condițiile atmosferice: straturile inferioare ale atmosferei provoacă o întârziere a semnalului electromagnetic emis de sateliți; întârzierea depinde de poziția satelitului pe boltă (cu cât satelitul este mai aproape de orizont cu atât semnalul va străbate o distanță mai mare prin straturile inferioare ale atmosferei – deci întârzierea va crește).

Eroarea de ceas a receptorului: principiul de funcționare a poziționării prin GNSS presupune ca ceasul satelitului care emite semnalul și cel al receptorului să fie perfect sincronizate – diferențele de sincronizare vor afecta astfel acuratețea de poziționare

Efectul acestor factori este redus semnificativ dacă se folosește metoda de poziționare denumită DGPS. În acest caz, echipamentul GNSS este în legătură continuă cu un număr de stații de referință (puncte a căror poziție este cunoscută pe care sunt montate echipamente GNSS). Fiecare stație își poate determina eroarea de poziționare curentă, ca diferența dintre poziția calculată de echipamentul GNSS și cea cunoscută. Această diferență este transmisă ca și corecție receptoarelor (în acest caz, echipamentului GNSS montat pe platforma de măsurare prin LiDAR) care o folosesc pentru a îmbunătăți precizia de poziționare.

Folosind această metodă, eroarea de poziționare va depinde de doi factori:

numărul de stații de referință la care este conectat receptorul

distanța dintre receptor și stațiile de referință.

Detalii suplimentare referitoare erorile specifice tehnologiei de poziționare GNSS, inclusiv pentru metoda DGPS, se regăsesc în Konecny (2014). Autorul remarcă faptul că folosirea acestei tehnologii poate îmbunătăți acuratețea de poziționare a echipamentelor GNSS mobile până la nivelul a 10-15 centimetri.

Erori asociate componentei IMU

Echipamentul IMU (Inertial Measurement Unit) determină orientarea în spațiu a platformei de măsurare prin calcularea înclinării pe cele trei axe (longitudinală, transversală și de rotație). Erorile de determinare a acestor unghiuri au efect asupra preciziei de determinare a coordonatelor punctelor înregistrate. Spre exemplu, eroarea unghiulară pentru axa transversală are efect în principal asupra determinării coordonatei Y a punctelor (de-a lungul direcției de scanare). În mod similar, eroarea unghiulară longitudinală are efect asupra coordonatei X (de-a lungul direcției de zbor).

Efectul erorilor unghiulare datorate componentei IMU crește odată cu unghiul de scanare și cu altitudinea de zbor (Baltsavias, 1999; Csanyi, 2007). În general, magnitudinea erorilor componentei IMU este estimată de producătorul echipamentului și inclusă ca parametru cunoscut în modelul de ajustare a datelor (Wang et al, 2008).

Calibrarea între senzori

Cele trei componente principale ale sistemului (senzor laser, IMU, GNSS) se află într-o poziție relativă fixă, care trebuie determinată cu precizie. În acest caz, apar două tipuri de erori:

Eroarea de măsurare a distanței relative dintre componente

Eroarea unghiulară dintre carcasa IMU și carcasa senzorului laser (eroare denumită și „nealinierea axelor”)

Prima eroare are valori de ordinul milimetrilor și un efect neglijabil asupra acurateții de poziționare a punctelor. În schimb, eroarea unghiulară este amplificată de altitudinea de zbor, astfel că efectul acesteia asupra acurateții înregistrărilor poate fi semnificativ (Csanyi, 2007).

Decalajul de timp

După cum s-a descris anterior, pentru poziționarea cu precizie a punctelor scanate este necesară determinarea următoarelor elemente: coordonatele platformei, unghiurile de rotație ale acesteia și distanța până la punctul aflat la sol. Pentru fiecare punct înregistrat, aceste elemente trebuie măsurate în același timp. Atunci când există un decalaj de timp între cele trei componente care le măsoară (emițător laser, IMU și GNSS), acesta trebuie determinat cu precizie și luat în calcul la procesarea datelor.

Diferența de timp atribuită componentei IMU este mai puțin importantă, întrucât unghiurile de rotație ale platformei sunt relativ stabile în timp. Totuși, această eroare influențează acuratețea de măsurare în condițiile unui zbor turbulent (Baltsavias, 1999).

Divergența razelor laser

Pentru simplificare, pulsul laser se modelează sub forma unei linii de grosime nulă. În realitate, asemenea oricărei unde electromagnetice care se propagă în aer, pulsurile laser sunt afectate de fenomenul de divergență. Astfel, presupunând cazul simplificat al unei suprafețe de teren orizontale, proiecția la sol a pulsului laser nu va fi un punct, ci un cerc (intersecția dintre conul rezultat prin divergența undelor electromagnetice și suprafața terenului). Diametrul acestei proiecții (denumită uzual „amprentă”) se poate determina cu formula:

[4.1]

unde:

D este diametrul amprentei pulsului laser;

r distanța dintre emițător și sol;

Δθ unghiul de divergență al razei;

Formula 4.1 se aplică pentru cazul simplificat al unei unde electromagnetice orientată pe verticală. Ca o corecție suplimentară, se poate modela efectul unei raze laser non-verticale precum și cel al terenului înclinat. În acest caz, diametrul amprentei se poate calcula cu formula:

[4.2]

unde:

β este unghiul de scanare al razei laser față de verticală;

α unghiul dintre normala suprafeței și verticala locului;

Mărimea amprentei la sol poate introduce erori semnificative de poziționare verticală ale punctelor, în special în cazul terenurilor înclinate (figura 4.4). Unghiul de divergență al razelor este o caracteristică de construcție a echipamentului de măsurare; astfel controlarea dimensiunii amprentei la sol se poate face prin reducerea altitudinii de zbor și limitarea scanării la unghiuri apropiate de verticală (cel mult 20-30 de grade).

În literatura de specialitate nu există un consens asupra valorilor acceptate pentru tipurile principale de erori, dată fiind evoluția constantă a tehnologiei. Pentru referință, în Tabelul 4.1 se prezintă valorile tipice pentru o parte dintre erori, așa cum au fost comunicate de către Shan și Toth (2009).

Fig. 4.4. Efectul amprentei la sol a razei laser în cazul terenurilor înclinate. Eroarea de determinare a distanței Δr va avea efect asupra acurateții de poziționare verticală a punctului.

Tabelul 4.1. Principalele surse de erori care afectează acuratețea măsurătorilor prin tehnologie ALS (după Shan și Toth, 2009).

4.3. Procesarea datelor colectate prin scanare laser aeriană

Transformarea datelor brute înregistrate de echipamentul pentru scanarea laser într-un set de date final, care poate fi utilizat în practică, implică următoarele etape principale:

Ajustarea între benzi

Deplasarea pe două axe a echipamentului de scanare (mișcarea longitudinală asigurată de deplasarea aeropurtată, iar mișcarea senzorului în plan transversal dată de mecanismul de baleiere) produce un tipar de înregistrare a punctelor sub forma unor benzi de scanare. Acesta depinde de următorii parametri:

Unghiul maxim de scanare.

Frecvența de emise a pulsurilor laser.

Viteza și altitudinea de zbor.

O anumită suprafață de teren va fi scanată de două sau mai multe ori, fiind cuprinsă în cel puțin două benzi învecinate. Întrucât măsurătorile sunt afectate de o serie de erori sistematice și întâmplătoare (descrise anterior), vor rezulta discrepanțe între benzile învecinate. Un anumit element din teren poate fi poziționat diferit, în funcție de banda din care acesta a fost scanat.

Aceste diferențe trebuie reduse la minim, pentru a asigura calitatea produsului final. Faptul că discrepanțele urmează un tipar sistematic permite modelarea și deci corectarea acestora. Procesul prin care se realizează această reducere a diferențelor între benzi se numește ajustare.

Până în prezent s-au dezvoltat mai mulți algoritmi de ajustare între benzi. Indiferent de algoritm, se urmărește determinarea unor elemente comune între benzi parțial suprapuse și calcularea discrepanțelor între acestea (erorile de poziționare în plan și verticale).

În domeniul fotogrammetriei, elementele comune pe baza cărora se determină diferențele dintre imaginile suprapuse sunt relativ ușor de identificat. În cazul scanării laser situația este îngreunată, deoarece punctelor nu li se asociază informație semantică care să permită identificarea lor fizică (Schenck, 2001). Unui punct scanat i se asociază doar date privind poziția și reflectanța, nu și o semnificație precum „cap de pod”, „colț clădire” sau „arbore”. Astfel este aproape imposibilă identificarea punctelor conjugate (Lee ș.a., 2007). Din acest motiv, algoritmii de ajustare între benzi urmăresc de regulă identificarea unor linii sau suprafețe comune (Habib ș.a., 2007).

Modul în care se folosesc aceste diferențe între elemente comune împarte algoritmii în două categorii (Shan și Toth, 2009):

Prima categorie urmărește modificări de tip „rubbersheet” aplicate asupra benzilor, pentru a reduce la minim discrepanțele dintre acestea; acești algoritmi nu sunt bazați pe modele fizice, urmărind strict reducerea erorilor de poziționare a punctelor. Algoritmii inițial dezvoltați au fost unidimensionali, urmărind doar ajustarea pe verticală (Kager și Kraus, 2001; Kornus și Ruiz, 2003). La scurt timp s-a descoperit importanța ajustării pe trei dimensiuni, adică diminuarea erorilor atât în plan vertical cât și orizontal (Maas, 2002; Vosselman, 2002b).

A doua categorie de algoritmi urmărește identificarea sursei erorilor, reducerea acestora realizându-se prin modificarea parametrilor senzorului și ale sistemului inerțial de orientare (practic se face o calibrare in-situ a echipamentului). Primul model de acest fel a fost propus de Behan et al (2000), cu alți algoritmi dezvoltați ulterior de Burman (2002), Filin (2003) sau Skaloud și Lichti (2006).

Clasificarea norului de puncte

Sistemul de măsurare LiDAR realizează scanarea suprafețelor prin efectuarea de măsurători stocate sub formă de puncte. Rezultă în urma zborului o mulțime de puncte (de ordinul milioanelor) poziționate într-un sistem de referință global, cunoscută uzual sub denumirea de nor de puncte.

Clasificarea este procesul prin care fiecărui punct i se atribuie ca etichetă o anumită categorie, în funcție de semnificația punctului respectiv (teren, clădire, vegetație etc.). După cum s-a descris anterior, un anumit puls laser are capacitatea de a produce una sau mai multe reflexii (sau „întoarceri”) pe traiectoria sa spre suprafața terestră. Atunci când deasupra suprafeței terestre se află anumite obstacole, acestea produc reflexii distincte. În cazul anumitor obstacole (precum vegetația, stâlpi sau fire de înaltă tensiune) doar o parte din puls este reflectat, astfel că ulterior se va produce o a doua reflexie. Alte obstacole (spre exemplu clădirile) reflectă întreaga energie a pulsului laser, împiedicând reflexii ulterioare.

Norul de puncte rezultat în urma ajustării între benzi va cuprinde astfel reflexii produse de diferite elemente din teren (figura 4.5).

Fig. 4.5. Secțiune printr-un nor de puncte înregistrat prin scanare laser aeriană.

Se observă reflexiile produse de suprafața terenului, dar și diferite obstacole aflate deasupra acesteia. Cu nuanțe de gri se reprezintă intensitatea semnalului reflectat.

De cele mai multe ori, în funcție de scopul colectării datelor, este necesară determinarea tipului de obiect asociat fiecărui punct (figura 4.6). Acest proces de clasificare se poate face:

manual, de către operatori specializați care analizează vizual norul de puncte împărțit pe profile și atribuie fiecărui punct categoria probabilă;

automat, caz în care se folosesc algoritmi de clasificare care iau în calcul informații precum: intensitatea reflexiei, altitudinea punctelor sau structură geometrică a norului de puncte;

combinat: rezultatul clasificării automate este analiza de operatori și eventualele erori identificate sunt corectate manual;

Stocarea datelor obținute prin scanare laser aeriană se face de regulă folosind formatul de fișiere .las, introdus de ASPRS. Acesta permite stocarea a maxim 655 clase distincte pentru puncte. 15 dintre aceste clase, cele mai uzuale, sunt prezentate în tabelul 4.2.

Fig. 4.6. Aceeași secțiune ca în figura 4.5, după clasificarea punctelor. Se observă suprafața terenului (maro), vegetație joasă (verde închis) sau medie (verde deschis), stâlp și cabluri de tensiune (albastru, respectiv galben) precum și o clădire (portocaliu).

Tabelul 4.2. Clasele de obiecte ce pot fi asociate punctelor LiDAR stocate în format .las (sursa: www.asprs.org).

Se fac următoarele mențiuni:

Valorile de clasă 0, respectiv 1 nu reprezintă obiecte în sine. Acestea identifică punctele care nu au fost clasificate (valoarea 0), respectiv punctele pentru care nu s-a putut identifica natura obiectului (valoarea 1 – neatribuit)

Valorile lipsă (8, 12) precum și valorile din intervalul 19-63 sunt rezervate de ASPRS

Valorile din intervalul 64-255 pot fi folosite pentru clase adiționale, stabilite de utilizatori

Clasa 7 (punct jos) cuprinde puncte cu o altitudine foarte joasă relativ la punctele învecinate, ceea ce înseamnă că sunt probabil înregistrări eroante. Similar pentru punctele din clasa 18 (zgomot), care au o altitudine relativ mare față de punctele învecinate. Ulterior clasificării, punctele asociate acestor două clase sunt de regulă eliminate din setul de date.

5. FILTRAREA DATELOR LIDAR

5.1. Aspecte generale

Una dintre etapele importante de prelucrare a datelor LiDAR este clasificarea norului de puncte. Atunci când scopul final este modelarea suprafeței terestre, avem de-a face cu un caz particular de clasificare, de tip binar: separarea punctelor aflate la sol de cele aflate deasupra terenului. Acest tip de clasificare poartă numele de filtrare a terenului, sau mai simplu – filtrare (Sithole și Vosselman, 2005; Chen et al, 2007).

Scopul urmărit prin filtrarea datelor este îndepărtarea acelor observații din norul de puncte inițial care se află deasupra terenului. Aceasta deoarece prezența acelor puncte nu ar permite modelarea corectă a suprafeței terestre. Astfel procesul de filtrare înfluențează în mod direct acuratețea de modelare a terenului.

În cazul terenurilor acoperite cu vegetație forestieră, proporția pulsurilor laser care penetrează întregul coronament, reușind să atingă suprafața terestră, este redusă. Din acest motiv, majoritatea punctelor înregistrate se vor afla deasupra terenului. Astfel, procesul de filtrare a datelor LiDAR brute este esențial pentru aplicațiile din domeniul forestier (Nord-Larsen și Riis-Nielsen, 2010; Maltamo et al, 2014).

Suprafața terestră, privită la nivel de detaliu, are de regulă un grad ridicat de complexitate. Astfel, pentru a se permite automatizarea filtrării datelor LiDAR, este necesară adoptarea unor presupuneri care simplifică realitatea (după Meng et al, 2010):

În cadrul unei anumite vecinătăți, punctul cu cea mai joasă elevație se află la sol – pornind de la această premisă, anumiți algoritmi de filtrare împart zona de lucru într-o rețea, folosind punctele de elevație minimă din fiecare celulă ca inițializare a filtrării (Masaharu și Ohtsubo, 2002; Silvan-Cardenas și Wang, 2006; Meng et al, 2009).

Panta între un punct aflat la sol și unul aflat deasupra solului este mai mare decât panta dintre două puncte la sol – astfel, o clasă de algoritmi bazați pe pantă determină panta între perechi de puncte – atunci când aceasta depășește o valoare limită, punctul cu altitudinea mai mare este eliminat (Vosselman, 2000; Sithole, 2001; Meng, 2005; Filin și Pfeifer, 2006).

Suprafața terenului este relativ netedă, fără discontinuități majore.

Se remarcă faptul că regulile detaliate mai sus nu sunt valabile în mod general. De exemplu, cea de-a doua presupunere își pierde din utilitate în cazul terenurilor accidentate.

5.2. Cercetări privind filtrarea datelor LiDAR

De-a lungul timpului s-au dezvoltat numeroși algoritmi de filtrare a datelor LiDAR, pentru domenii de activitate sau condiții de teren diverse; totodată, se dezvoltă metode de comparare a performanțelor acestor algoritmi, pentru situații cât mai variate.

Prima modalitate automată, algoritmică de filtrare a datelor obținute prin măsurători LiDAR este propusă de Kraus (1997). În anul următor, Kraus și Pfeifer (1998) detaliază o metodologie de interpolare/filtrare a măsurătorilor LiDAR, adresată modelării suprafețelor acoperite de vegetație forestieră. Axelsson (2000) propune un algoritm de filtrare bazat pe generarea de Rețele Neregulate de Triunghiuri (TIN). Metoda presupune generarea unui TIN pe baza unor puncte de inițializare, rețea care reprezintă o primă aproximare grosieră a suprafeței terestre. Punctele sunt clasificate în funcție de distanța și unghiul față de cea mai apropiată fațetă a rețelei. În cazul în care cele două elemente geometrice depășesc valori-prag stabilite, punctul este eliminat. În caz contrar, se consideră că punctul este localizat la suprafața terenului, fiind inclus într-un TIN ulterior. Astfel, rețelele generate iterativ se apropie succesiv de suprafața reală a terenului. În prezent clasa de algoritmi care urmează aceste principii (cunoscută sub numele de PTD) se bucură de popularitate, pentru acuratețe, eficiență dar și datorită integrării în programul comercial de procesare LiDAR TerrascanTM al companiei Terrasolid (Zhao et al, 2016).

Terenurile împădurite reprezintă o situație aparte pentru filtrare, dată find densitatea punctelor la sol semnificativ redusă (prin comparație cu terenurile descoperite) în combinație cu topografia complexă care se asociază în general cu vegetația forestieră (Maguya, 2014). Din acest motiv, una dintre direcțiile de cercetare este dezvoltarea algoritmilor de filtrare pentru terenurile împădurite. Dintre acești algoritmi amintim: Groundfilter (Kraus și Pfeifer, 1998), MCC (Evans și Hudak, 2007), REIN (Kobler et al, 2007), IPTD (Zhao et al, 2016) sau algoritmii dezvoltați de Maguya et al (2014), Guan et al (2014), Kim și Eo (2015).

Dintre direcțiile de cercetare recente, amintim dezvoltarea de algoritmi morfologici (Zhang et al, 2003; Arefi și Hahn, 2005; Chen et al, 2007; Pingel, 2013), algoritmi bazați pe segmentarea norului de puncte (Tovari și Pfeifer, 2005; Zhang și Lin, 2013) și algoritmi care folosesc învățarea automată sau rețele neurale (Hu și Yan, 2016).

Odată cu înmulțirea numărului de algoritmi de filtrare apare și necesitatea comparării acestora din punct de vedere al acurateții. Sithole și Vosselman (2004) delimitează opt zone-test, din mediul urban și rural, care acoperă o gamă largă de folosințe ale terenului (suprafețe acoperite de pădure, clădiri, albii traversate de poduri etc.). Autorii realizează o comparație calitativă și cantitativă a opt algoritmi de filtrare pentru aceste suprafețe, identificând condițiile dificile pentru filtrare. Zhang și Whitman (2005) compară performanțele a trei algoritmi pentru filtrarea a trei suprafețe-test, una dintre ele fiind caracterizată ca “pădure cu relief accidentat”. Korzeniowska et al (2014) testează cinci algoritmi, în condiții variate de topografie și folosință a terenului. Se remarcă procedura formală de analiză a acurateții, acoperirea spațială a diferitelor tipuri de erori fiind estimată și comparată. Montealegre et al (2015) testează șapte algoritmi de tip sursă-liberă în condițiile unei păduri de tip mediteranean. Autorii investighează influența unor factori precum panta terenului, categoria de folosință sau densitatea norului de puncte.

Este important de menționat că numeroasele studii comparative ale acurateții algoritmilor de filtrare nu au condus la un consens. Variabilitatea condițiilor de teren dar și a datelor inițiale (densitatea punctelor, tipul de senzor, altitudinea de zbor etc.) îngreunează formularea unor concluzii generale. Totuși, studiile realizate până în prezent identifică o serie de condiții care pun probleme tuturor algoritmilor de filtrare, spre care este deci de interes orientarea cercetărilor viitoare. Acestea sunt: relief accidentat, discontinuități ale pantei sau prezența vegetației forestiere dense, în special pe terenuri cu înclinație ridicată (Meng 2010; Pingel, 2013).

5.3. Tipuri de erori la filtrare

Prin filtrarea norului de puncte înregistrate folosind tehnologia LiDAR se obține un set de puncte redus. Acest subset al norului de puncte inițial include acele observații care au fost considerate ca reprezentând suprafața terestră.

În cadrul acestui subset de observații se întâlnesc două tipuri de erori, care influențează acuratețea de modelare ulterioară a suprafeței terestre:

a. Eroarea de comitere:

O observație inclusă în norul de puncte LiDAR, reprezentând un obstacol aflat deasupra terenului, poate să nu fie eliminată prin procedura de filtrare. Această observație poartă denumirea de eroare de comitere.

În acest caz, eroarea de determinare a altitudinii este una pozitivă (modelul suprafeței terestre va avea o altitudine mai mare decât cea reală). În figura 5.1 (a) se prezintă efectul unor astfel de erori asupra modelului suprafeței. Eroarea de comitere este cunoscută și sub numele de eroare de tip I.

b. Eroarea de omitere:

În cazul invers erorii de comitere, prin filtrare se poate elimina o observație din norul de puncte, deși în realitate aceasta se află la nivelul terenului. În acest caz, se poate vorbi de o eroare de omitere. Efectul acesteia asupra determinării altitudinii va fi unul negativ (modelul suprafeței terestre va avea o altitudine mai mică decât cea reală). În figura 5.1 (b) se prezintă efectul unor astfel de erori asupra modelului suprafeței. Eroarea de comitere este cunoscută și sub numele de eroare de tip II.

(a) (b)

Fig. 5.1. Efectul erorilor de filtrare de tip I (a), respectiv al erorilor de filtrare de tip II (b) asupra modelării terenului.

Interpolarea suprafeței, care are ca efect scăderea variației valorilor de altitudine poate produce o reducere a efectului erorilor de filtrare, îmbunătățind calitatea reprezentării finale a suprafeței terestre (Hodgson și Bresnahan, 2004).

6. INTERPOLAREA SUPRAFEȚEI TERESTRE PE BAZA MĂSURĂTORILOR LIDAR

6.1. Aspecte generale

Variabilitatea în spațiu a fenomenelor fizico-geografice se estimează prin eșantionarea unui număr finit de puncte, de regulă iregular distribuite în timp și spațiu. În același timp, vizualizarea, analiza sau modelarea acestor fenomene în cadrul Sistemelor de Informații Geografice se face prin reprezentarea sub formă de suprafețe continue. Transformarea unui set de puncte într-o suprafață continuă necesită estimarea valorilor pentru locații noi, ne-eșantionate (figura 6.1). Acest proces se numește interpolare sau aproximare. Interpolarea presupune determinarea valorii unui punct ne-eșantionat, pe baza măsurătorilor realizate pentru puncte învecinate acestuia, denumite puncte de referință (El-Sheimy ș.a., 2005). Prin contrast, extrapolarea este procesul de estimare pentru puncte situate în afara zonei de eșantionare.

În ceea ce privește suprafața terestră, o reprezentarea fidelă a acesteia ar necesita în principiu determinarea altitudinii pentru un număr infinit de puncte, imposibilă în practică (Tamaș și Tereșneu, 2010). În fapt, indiferent de metoda de înregistrare a datelor altimetrice, se fac determinări pentru un număr limitat de puncte, care stau la baza generării suprafețelor. Modelul suprafeței terestre este de regulă de tip raster, adică o rețea de celule egale ca dimensiune cărora li se atribuie o valoare singulară pentru altitudine. Acest model este cunoscut sub numele de Model Digital Altimetric (MDA). Generarea unui MDA pe baza punctelor eșantionate necesită estimarea, prin interpolare, a valorii coordonatei Z pentru centrele celulelor modelului, pe baza valorilor Z ale punctelor cunoscute din vecinătate.

Fig.6.1. Interpolarea valorii pentru punctul A

pe baza valorilor cunoscute pentru punctele 1-4.

6.2. Generarea suprafețelor prin interpolare

Procesul de interpolare este unul matematic, deci se poate realiza automat, folosind algoritmi specifici. Aceștia urmăresc definirea unei funcții de predicție, cu ajutorul căreia se estimează valorile pentru puncte noi (care nu sunt incluse în setul de puncte eșantionate). Forma generală a acestei funcții este (Webster și Oliver, 2001):

[6.1]

unde:

este valoarea estimată pentru punctul x0;

z valoarea cunoscută pentru punctul xi;

λi ponderea atribuită observației z;

n numărul de puncte cunoscute folosit la estimarea valorii lui ;

Un algoritm de interpolare poate fi (Li și Heap, 2014):

determinist (sau non-geostatistic): parametrii modelului sunt determinați în mod empiric sau arbitrar și nu se obține o estimare a erorilor de modelare;

probabilistic (sau geostatistic): parametrii modelului sunt determinați în mod obiectiv, iar în urma modelării rezultă o estimare a distribuției spațiale a erorilor de modelare;

De asemenea, algoritmii se pot clasifica în următoarele categorii (după Hengl, 2009):

exacți – pentru un punct cunoscut, valoarea prezisă este egală cu cea măsurată (cu alte cuvinte, suprafața interpolată va trece prin toate punctele de referință); aproximativi -valoarea prezisă nu este în mod necesar egală cu cea măsurată;

locali – doar punctele cunoscute din proximitatea locației au influență asupra predicției; globali – toate punctele măsurate influențează predicțiile;

convecși – valorile prezise se încadrează în intervalul valorilor cunoscute; non-convecși valorile prezise pot depăși intervalul valorilor cunoscute;

Detalii suplimentare privind metodele de interpolare folosite în geomatică se regăsesc în Wood și Fisher (1993), MacMillan et. al (2008), Li și Heap (2008) sau Lam (2013).

6.3. Erori de interpolare a măsurătorilor LiDAR

Interpolarea este procesul prin care se atribuie valori pentru un set de puncte noi, pe baza valorilor punctelor eșantionate învecinate. Procesul este unul aproximativ, prin care valorile sunt estimate, fiind deci afectat de erori de determinare. O parte din cercetările efectuate asupra erorilor de interpolare a măsurătorilor LiDAR sunt rezumate în continuare.

Hodgson și Bresnahan (2004) acordă erorii de interpolare a doua poziție ca importanță (după eroarea sistemului de măsurare), într-o clasificare a componentelor erorii totale a Modelelor Altimetrice Digitale derivate din măsurători LiDAR. Numeroase studii empirice au urmărit estimarea efectului incertitudinii de interpolare asupra acurateții produsului final. În continuare vor fi descrise sumar o parte dintre acestea, toate referitoare la analiza erorilor de interpolare pentru cazul specific al datelor de altimetrie obținute prin măsurători LiDAR.

Hyyppa et al (2005) testează un algoritm de interpolare propriu în condițiile unei zone de studiu împădurite. Autorii estimează o eroare standard de 0.15 m pentru terenurile relativ plane, respectiv o eroare standard de 0.40 m pentru terenurile cu pante de peste 40 de procente.

Su și Bork (2006) compară algoritmii de interpolare Inverse Distance Weighted (IDW), Kriging și spline în condițiile unui relief complex (pe care îl denumesc „vălurit”) din cadrul Aspen Parkland, Canada. Autorii concluzionează că erorile de interpolare sunt cele mai mici în cazul algoritmului IDW, urmat de Kriging. Acuratețea modelului este afectată de pantă, eroarea medie pătratică (RMSE) fiind dublă în zonele cu pante de peste 15 grade, comparativ cu zonele cu pante de sub 2 grade. De asemenea, autorii constată că erorile depind și de acoperirea cu vegetație: pentru terenurile împădurite elevația este supra-estimată iar pentru pajisți elevația este sub-estimată.

Bater și Coops (2009) evaluează erorile de interpolare pentru algoritmii Inverse Distance Weighted (IDW), Natural Neighbour, Linear, Spline, Quintic și ANUDEM. Se analizează ce efect au asupra erorilor următorii factori: rezoluția modelului (rezoluțiile testate: 0.5, 1.0 și 1.5 metri), densitatea punctelor, înclinarea terenului și clasa de vegetație. Zona de studiu este predominant forestieră, cu păduri mature. Cu excepția algoritmului IDW, care oferă rezultate comparativ mai slabe, suprafețele interpolate de ceilalți algoritmi au acuratețe similară, atât global cât și pentru stratificarea în funcție de pantă sau densitate a punctelor. Odată ce parametrii algoritmilor sunt optimizați, autorii observă că rezoluția modelului afectează acuratețea, creșterea acesteia conducând la erori mai mari.

Ismail (2016) folosește puncte determinate cu stație totală și nivelă optică pentru a testa acuratețea MDA pentru suprafețe acoperite cu arbori de cauciuc și palmier pentru ulei. Algoritmii de interpolare testați sunt: Inverse Distance Weighted (IDW), Kriging și Spline. Rezultatele arată că performanța algoritmilor depinde de tipul de vegetație: pentru suprafața cu arbori de cauciuc algoritmul Spline are cea mai mică valoare pentru eroarea medie pătratică, în timp ce pentru suprafața acoperită cu palmier algoritmul Kriging este cel mai precis. Pentru fiecare algoritm erorile cresc odată cu creșterea pantei terenului și cu modificarea rezoluției modelului (de la 0.5 la 1.0 metri).

Sterenczak et al (2016) înregistrează date de altimetrie folosind stație totală și măsurători GPS RTK (Real Time Kinematic) pentru a testa acuratețea de interpolare în cazul suprafețelor montane împădurite. Autorii testează interpolarea cu Inverse Distance Weighted (IDW), Triangular Irregular Network (TIN) și Natural Neighbour. Se testează de asemenea o serie de algoritmi implementați în programe specializate de procesare a datelor LiDAR (TerraScanTM, Opals și TreesVis). Cele mai precise suprafețe rezultă prin folosirea algoritmilor Natural Neighbour sau TIN, însă autorii conchid că toate metodele de interpolare conduc la erori similare dacă se aleg parametrii optimi pentru procesare. Înclinarea terenului și prezența vegetației tinere sunt identificați ca factorii cu influență asupra preciziei de modelare a suprafeței terestre.

7. ISTORIC AL CERCETĂRILOR PRIVIND SCANAREA LASER ÎN DOMENIUL FORESTIER

Dezvoltarea tehnologiei de măsurare care folosește impulsuri electromagnetice de tip laser a avut loc concomitent cu realizarea unor progrese importante în domenii variate, precum:

dezvoltarea și răspândirea Sistemelor de Informații Geografice;

implementarea primului sistem de poziționare global de către SUA;

dezvoltarea tehnologiei de poziționare folosind echipamente inerțiale (IMU – Inertial Measurement Unit);

apariția calculatoarelor de tip PC (computer personal) și evoluția capacităților de stocare și procesare a datelor;

Acest cumul de factori a permis răspândirea sistemelor de măsurători laser în toate domeniile în care este necesară cartarea suprafeței terestre.

În silvicultură, potențialul tehnologiei este remarcat încă din anii 1990. Friedrich Ackermann, „părintele” răspândirii tehnologiei LiDAR în comunitatea academică europeană (Carson ș.a., 2004) remarcă faptul că motivația sa inițială pentru cercetarea în domeniul LiDAR a fost dată de nevoile domeniului forestier (ISPRS, 1999). De altfel, primele cercetări cu privire la automatizarea filtrării datelor LiDAR au vizat terenurile acoperite cu vegetație forestieră (Kraus și Pfeifer, 1998). În prezent se poate spune că tehnologia LiDAR a adus o nouă percepție asupra structurii tridimensionale a pădurii, la nivel local sau regional (Bauwens et al, 2016).

Potențialul scanării laser pentru domeniul forestier este vast, tehnologia permițând măsurarea unor caracteristici importante ale pădurii. Înregistrarea punctelor la nivelul superior al coronamentului, la nivelul suprafeței terestre, precum și între cele două nivele permite modelarea structurii verticale a pădurii. Primele cercetări privind estimarea structurii arboretelor folosind date LiDAR au fost realizate în 1986 (Maclean și Krabill, 1986 – citat în Lefsky et al, 2002). Cercetări ulterioare au dezvoltat în continuare domeniul, prin sporirea numărului de parametri ce pot fi estimați, îmbunătățirea acurateții de modelare a algoritmilor sau testarea în diferite ecosisteme (Popescu et al, 2002; Clark et al, 2004; Goodwin et al, 2006; Van Leeuwen et al, 2010; Popescu et al, 2011; Birjaru, 2011; Ferraz et al, 2016).

Dintre cele mai importante ramuri de cercetare sau aplicative ale silviculturii în care tehnologia LiDAR și-a adus aportul, amintim:

Inventarierea pădurilor: primele cercetări privind folosirea tehnologiei LiDAR pentru inventarierea pădurilor au avut loc cu precădere în țările scandinave, începând din anii 1990 (Naesset, 1997; Hyyppa ș.a., 2009); ulterior tehnologia s-a răspândit în restul Europei și pe continentul american (Wallace ș.a., 2012).

Monitorizarea stării de sănătate a pădurii: deși imaginile satelitare oferite de programele Landsat, MODIS sau SPOT rămân încă cea mai populară sursă de date pentru monitorizarea extensivă a pădurii, Rullan-Silva ș.a. (2013) remarcă faptul că tehnologia LiDAR devine tot mai populară în acest domeniu (Solberg ș.a., 2004; Vastaranta ș.a., 2013).

Identificarea arborilor individuali: densitatea ridicată a observațiilor LiDAR permite în prezent cartografierea structurii pădurii la nivel de arbori individuali, fără necesitatea deplasării în teren (Popescu și Wynne, 2004; Chen et al, 2006; Jakubowski et al, 2013; Wu et al, 2016).

Estimarea cantității de biomasă: în acest domeniu de importanță tot mai mare (datorită schimbărilor climatice recente) măsurătorile LiDAR au o tradiție îndelungată, fiind folosite încă din anii 1980 (Nelson et al, 1988; Patenaude et al, 2004; Popescu et al, 2011; Ferraz et al, 2016).

În România, tehnologia LiDAR poate fi considerată încă în faza de pionerat. Dintre cercetările realizate în silvicultură, remarcăm folosirea scanării laser pentru: estimarea caracteristicilor structurale ale pădurii (Birjaru, 2011; Petrila ș.a., 2012a; Apostol, 2016), estimarea volumului de masă lemnoasă pe picior pentru arborete de molid (Apostol ș.a., 2012), estimarea cantității de biomasă (Petrila ș.a., 2012b) sau detecția coroanelor (Strâmbu și Strâmbu, 2015).

Majoritatea parametrilor structurali estimați pe baza măsurătorilor LiDAR (de ex. înălțimea coronamentului, înălțimea arborilor individuali, cantitatea de biomasă) necesită modelarea a priori a suprafeței terestre (Tinkham et al, 2012). Astfel, una din direcțiile de cercetare din domeniu este dezvoltarea algoritmilor de filtrare a măsurătorilor LiDAR, pentru cazul specific al terenurilor acoperite cu vegetație forestieră. Descrierea principalilor algoritmi de filtrare din această categorie se regăsește în secțiunea 9.4.4.

Analiza surselor de erori de altimetrie specifice măsurătorilor LiDAR constituie de asemenea o direcție de cercetare importantă (Hopkins et al, 2004; Hodgson et al, 2005; Aguilar et al, 2010; Tinkham et al, 2012). Dintre studiile care au urmărit evaluarea acurateții de reprezentare a suprafeței terestre în cazul terenurilor împădurite, menționăm următoarele: Hyppa ș.a. (2000), Reutebuch ș.a. (2003), Hodgson și Bresnahan (2004), Naesset ș.a. (2015), Sterenczak ș.a. (2016), Cățeanu și Arcadie (2017).

Pe lângă faptul că permite determinarea parametrilor structurali ai pădurii, reprezentarea suprafeței terestre poate servi în mod direct la aplicații din domeniul gospodăririi sau exploatării pădurilor. Dintre acestea, amintim posibilitatea identificării sau proiectării drumurilor forestiere (Aruga ș.a., 2005; Pentek ș.a., 2005; White ș.a., 2010), amenajarea căilor de colectare a lemnului (Sterenczak și Moskalik, 2014), alegerea soluțiilor de colectare (Mohtashami ș.a., 2012; Vega-Nieva ș.a., 2009) sau planificarea operațiunilor de exploatare (Chung, 2003).

În privința tendințelor de cercetare recente, se remarcă fuzionarea datelor LiDAR cu alte surse de date (de ex. imagini capturate cu echipamente de tip UAV) pentru o modelare mai precisă (Hyde et al, 2006; Dalponte et al; 2012;Aguilar et al, 2019), modelarea de la sol a structurii pădurii folosind echipamente de scanare fixe (Moskal și Zeng, 2012; Stovall ș.a., 2017; Schneider ș.a., 2019) sau echipamente instalate pe platforme mobile (Strahler ș.a., 2008; Tang ș.a., 2015; Bauwens ș.a., 2016; Liang ș.a., 2018).

8. CARACTERISTICI ALE CONDIȚIILOR DE CERCETARE

8.1. Scopul și obiectivele cercetărilor

Definirea obiectivelor cercetării a pornit de la posibilitatea de obținere a datelor LiDAR înregistrate prin tehnologia scanării laser aeriene (ALS), pentru suprafețe de teren împădurite și preferabil cu un grad ridicat de variabilitate sub aspect topografic.

În acest context, scopul cercetărilor este evaluarea potențialului tehnologiei LiDAR ca sursă de date pentru reprezentarea suprafeței terestre la nivel de detaliu, în domeniul silviculturii.

Pentru atingerea scopului propus, s-au stabilit o serie de obiective:

Cunoașterea tehnologiei de scanare cu echipament LiDAR, din perspectivă teoretică și ca aplicabilitate practică generală, dar mai ales în domeniul silviculturii.

Realizarea unei bibliografii extinse, cu referire la tehnicile geomaticii forestiere, algoritmii de procesare a datelor LiDAR, studii experimentale privind integrarea tehnologiei LiDAR în domeniul forestier ș.a.

Obținerea datelor LiDAR, procesarea preliminară a acestora și delimitarea zonelor de studiu.

Efectuarea de măsurători geo-topografice în teren, pentru obținerea datelor independente necesare evaluării preciziei.

Identificarea acelor etape de procesare a datelor LiDAR care influențează acuratețea de modelare a suprafeței terestre.

Procesarea datelor LiDAR, cu scopul final de a obține modele ale suprafeței terestre care pot analizate din punct de vedere al acurateții altimetrice.

Stabilirea unei metodologii de identificare, evaluare și analizare a erorilor produse de procesarea datelor, pe cât se poate obiectivă și care să permită analiza comparativă a metodelor de procesare studiat.

Interpretarea și analiza rezultatelor, din prisma distribuției statistice și spațiale a erorilor; estimarea efectului factorilor externi (condiții de teren, vegetație etc.) asupra variației erorilor.

Elaborarea concluziilor, urmată de sugestii de cercetare viitoare.

8.2. Locul cercetărilor

La noi în țară tehnologia LiDAR nu a cunoscut încă o largă răspândire în domeniul forestier, în principal datorită costurilor ridicate și necesarului de specialiști care să o poată pune în aplicare (Birjaru, 2011). După cunoștințele noastre, această afirmație rămâne valabilă și la momentul elaborării acestei lucrări.

În aceste condiții, singura soluție disponibilă pentru atingerea obiectivelor tezei a fost utilizarea unui set de date LiDAR înregistrate pentru Agenția Națională Apele Române, pus la dispoziție prin amabilitatea Institutului Național de Cercetare-Dezvoltare în Silvicultură „Marin Drăcea”.

Ținând seama de situația expusă anterior, alegerea zonei de studiu s-a făcut pe baza criteriului disponibilității datelor care pot fi utilizate pentru scopurile de cercetare propuse. Astfel, s-au stabilit două areale de cercetare (figura 8.1):

în Carpații Meridionali, în aproprierea localității Mălaia, județul Vâlcea;

în Carpații Meridionali, în apropierea Lacului de acumulare Vidra, județul Vâlcea;

În ambele locații au fost realizate măsurători LiDAR cu ajutorul platformelor aeropurtate, pentru gestionarea zonelor de risc limitrofe cursurilor de apă (în principal riscul de eroziune). Suprafețele pe care s-au făcut înregistrările cuprind zone împădurite, fiind caracterizate totodată de un relief accidentat, specific regiunii. Datele puse la dispoziție au fost colectate în 2008 și 2011 (pentru arealul Mălaia), respectiv în anul 2012 (pentru arealul Vidra).

Din punct de vedere geografic, arealul Mălaia se înscrie în unitatea de relief Carpații Meridionali, regiunea Munții Căpățânii. Forma de relief predominantă este versantul, iar panta terenului are valori între 20-30 de grade pe cea mai mare parte a suprafeței, atingând valori maxime de 46 de grade. Altitudinea minimă este de 480 metri, iar cea maximă de 1000 metri.

Sub aspectul vegetației forestiere, speciile principale sunt molidul și fagul, în arborete pure sau în amestecuri. Vârsta medie a arboretelor de fag este de 60-90 de ani, iar înălțimea medie este de 25-30 metri. Arboretele de molid au o vârstă medie între 80 și 100 de ani, iar înălțimea medie este în jur de 23-25 metri, conform datelor de amenajament.

În ceea ce privește suprafața Vidra, aceasta se încadrează în regiunea Munții Latoriței, parte a Carpaților Meridionali. Forma de relief predominantă este de asemenea versantul, terenul având înclinația între 5 și 35 de grade pe cea mai mare parte a suprafeței. Altitudinea are valori între 1300 și 1900 metri, conform datelor de amenajament.

Suprafața este încadrată în U.P. IV – Puru, Ocolul Silvic Voineasa, speciile principale fiind molidul și, într-o mai mică măsură, fagul. Zona sudică a suprafeței este acoperită cu pajiști alpine.

Fig. 8.1. Localizarea suprafețelor cu măsurători LiDAR disponibile

8.3. Baza de date

8.3.1. Date LiDAR

Datele LiDAR reprezintă suportul de bază în cercetarea de față. Cele două areale descrise anterior au fost cartate prin scanare laser aeriană (ALS – Airborne Laser Scanning), rezultând nori de puncte poziționați 3D care descriu suprafața terestră precum și obiectele aflate pe aceasta.

Datele au fost puse la dispoziție de către INCDS „Marin Drăcea” în format .las, în formă prelucrată (după ajustarea între benzi) și georeferențiate în sistemul UTM (Universal Transversal Mercator) zona 35N, datum S-42 Pulkovo.

În ceea ce privește clasificarea norului de pentru, pentru arealul Mălaia doar o parte din suprafață este clasificată, restul fiind oferită sub formă de fișiere .las fără clasificare. Pentru arealul Vidra, întreg norul de puncte pus la dispoziție este clasificat.

Analiza vizuală a datelor a evidențiat un număr redus de erori de clasificare semnificative, care au fost corectate manual.

Pentru scanarea celor două suprafețe s-a folosit un sistem compus din: senzor laser Riegl LMS-Q560, echipament de poziționare globală de tip GPS și sistem inerțial de orientare (IMU). Sistemul a fost instalat pe un avion Diamond Aircraft Industries, model DA42 MPP.

Fig. 8.2. Delimitarea arealului de studiu Mălaia, suprapusă pe imagine satelitară Bing Maps

(© Microsoft Corporation 2019).

Fig. 8.3 Delimitarea arealului de studiu Vidra, suprapusă pe imagine satelitară Bing Maps (© Microsoft Corporation 2019)

8.3.2. Date suport

Pe lângă datele LiDAR, pentru îndeplinirea obiectivelor propuse au fost folosite o serie de date auxiliare, respectiv:

planuri topografice având curbe de nivel (scara 1:25 000)

ortofotoplanuri

hărți gratuite din cadrul proiectului OSM (Open Street Map)

imagini satelitare de rezoluție medie, disponibile pe platformele Bing Maps (© Microsoft Corporation 2019) sau Google Maps (© Google 2019)

Model Digital de Elevație SRTM, disponibil gratuit pe platforma Earth Explorer a USGS (United States Geological Survey), la rezoluție de 1 arc-secundă (aprox. 38 metri în proiecție UTM)

Aceste date suport au servit la localizarea zonelor de interes, la delimitarea suprafețelor de probă, precum și la planificarea și desfășurarea etapei de teren.

8.4. Echipamente și instrumente pentru cercetare

8.4.1. Sistemul de scanare laser

Datele LiDAR au fost înregistrate folosind un sistem de scanare laser cu o structură convențională (senzor laser, sistem de poziționare globală de tip GPS și sistem de măsurători inerțiale – IMU) instalat la bordul unei aeronave Diamond Aircraft Industries.

Echipamentul laser model Riegl LMS-Q560 este un sistem cu măsurare continuă, echipat cu oglindă rotativă poligonală care produce un tipar de scanare sub forma unor linii paralele. Schema constructivă a sistemului se regăsește în figura 8.4, iar detalii tehnice adiționale sunt prezentate în tabelul 8.1.

Fig. 8.4. Componența sistemului de scanare laser LMS-Q560 dezvoltat de compania Riegl

(după broșură tehnică Riegl: www.riegl.com)

Tabelul 8.1. Parametrii senzorului LiDAR Riegl LMS-Q560 (sursă: LMS-Q560 Datasheet, Riegl).

8.4.2. Instrumente topografice și geodezice

Necesitatea realizării de măsurători în teren cu o precizie de poziționare comparabilă cu cea oferită de sistemul LiDAR a impus folosirea instrumentului de tip stație totală pentru cartarea terenului, respectiv a două echipamente de tip GNSS pentru geo-referențierea măsurătorilor. Instrumentele geotopografice au fost puse la dispoziție prin amabilitatea Facultății de Silvicultură și Exploatări Forestiere din cadrul Universității „Transilvania” Brașov.

Stația totală folosită, model Leica TCR 407 (figura 8.5), are o structură clasică, cu trei părți: componenta mecanică, optică și electronică. Dintre caracteristicile tehnice ale dispozitivului amintim:

acuratețe de determinare a unghiului: 7 secunde

acuratețe de determinare a distanțelor: 2mm + 2ppm (părți pe milion)

capacitatea de mărire a lentilei: 30x

distanță maximă de măsurare folosind unde laser: 3500 m

nivelă electronică

compensare electronică pe două axe, pentru reducerea erorilor de calare

autonomie de lucru sporită, datorită capacității de stocare de 10 000 înregistrări

Fig. 8.5. Stație totală Leica TCR 407, folosită la etapa măsurătorilor de teren.

Echipamentul de poziționare globală prin satelit, model Leica SR20, a fost folosit pentru georefențierea măsurătorilor realizate cu stația totală în sistemul de proiecție UTM 35N, permițând astfel suprapunerea măsurătorilor topografice cu datele LiDAR. Leica SR20 (figura 8.6.) este un echipament de tip GPS cu 12 canale, în simplă frecvență, care permite determinarea punctelor cu precizie centimetrică, pentru măsurători statice cu post-procesare.

Fig. 8.6. Echipament GPS Leica SR20, folosit la geo-referențierea măsurătorilor topografice.

8.4.3. Suportul informatic al cercetărilor

Pentru atingerea obiectivelor propuse, au fost necesare o serie de platforme software specializate, precum și generale. Suportul hardware al acestor platforme este reprezentat de un calculator personal, cu următoarele specificații:

unitate de procesare: AMD A10-6800K, 4 cores @ 4.10 Ghz

unitate de procesare grafică: Nvidia GTX 1050-Ti, 4 GB DDR5, 1392 MHz

memorie RAM: Kingston HyperX, 8 GB DDR3, 1333 Mhz

sistem de operare: Windows 10 Pro, arhitectură 64-bit; Linux KDE Neon arhitectură 64-bit

În ceea ce privește componenta software, se disting trei categorii principale de programe care au permis atingerea obiectivelor stabilite:

Programe pentru procesarea măsurătorilor topogeodezice:

Leica Geo Office, v. 8.3, pentru descărcarea datelor și post-procesarea măsurătorilor GPS;

Trimble Terramodel, v. 10.61, folosit la prelucrarea măsurătorilor topografice realizate cu stația totală;

Programe pentru procesarea și analiza datelor LiDAR:

pachetul software FUSION v. 3.50, dezvoltat de Robert J. McGaughey în cadrul Serviciului de Silvicultură al Departamentului de Agricultură al SUA (McGaughey, 2014); program disponibil în mod gratuit care permite procesarea, vizualizarea și analizarea datelor LiDAR pentru domeniul forestier;

suita de utilitare LAStools, dezvoltată de rapidlasso GmbH, care cuprinde numeroase unelte pentru procesarea, interpretarea, vizualizarea sau analizarea datelor LiDAR stocate în formatul .las;

Programe pentru Sisteme Informatice Geografice:

Global Mapper, v. 19, folosit la procesarea preliminară a datelor: schimbarea sistemului de proiecție, conversia între diverse formate de fișiere etc.;

GRASS GIS, v. 7.6.1, folosit în diferite etape de cercetare datorită posibilităților extinse de automatizare a operațiunilor oferite de limbajul de programare Python;

SAGA GIS, v. 2.3.2, cuprinde numeroase unelte pentru generarea, interpretarea și analizarea suprafețelor în format raster;

QGIS, v. 3.5, folosit în principal pentru vizualizarea datelor și generarea hărților;

Programe pentru analiza statistică:

Statistica, v. 10, folosit pentru analiza statistică a erorilor;

mediul software gratuit R, v. 3.5 (R Core Team, 2016), folosit la generarea graficelor;

9. METODOLOGIA DE CERCETARE

9.1. Prezentare generală

Unul dintre obiectivele cercetării a fost identificarea acelor etape de procesare a datelor LiDAR care influențează gradul de acuratețe al produsului final (Modelul Digital Altimetric). În urma analizării literaturii de specialitate s-au identificat trei factori principali cu efect asupra acurateții Modelului Digital Altimetric (MDA). Aceștia sunt:

Caracteristicile datelor de intrare: precizie de poziționare, densitatea punctelor, numărul de reflexii per puls etc.; acestea depind de elemente precum: senzorul laser, altitudinea de zbor, viteza de deplasare a platformei, unghiul maxim de scanare sau acoperirea terenului.

Procedura de filtrare a datelor LiDAR, prin care se urmărește separarea punctelor aflate la sol de punctele reprezentând obstacole aflate deasupra suprafeței terestre (în cazul de față obstacolul este de regulă vegetația forestieră).

Procedura de interpolare a setului de puncte filtrat, prin care se generează modelul suprafeței terestre.

Întrucât nu a fost posibilă realizarea unui campanii proprii de înregistrări LiDAR, nu s-a putut interveni asupra factorilor din prima categorie. Astfel, cercetarea s-a concentrat asupra ultimelor două categorii de factori care influențează acuratețea de modelare a suprafeței terestre.

Pentru atingerea obiectivelor propuse s-au definit următoarele etape de lucru:

Analiza preliminară a seturilor de date LiDAR disponibile și delimitarea unor suprafețe de studiu.

Prelucrarea preliminară a datelor LiDAR.

Planificarea și realizarea măsurătorilor în teren.

Definirea, calcularea și analiza erorilor produse de procesul de filtrare a datelor LiDAR.

Definirea, calcularea și analiza erorilor produse de procesul de interpolare a suprafeței.

9.2. Analiza preliminară a datelor LiDAR

După cum s-a descris anterior, pe baza datelor LiDAR avute la dispoziție, s-au stabilit două areale de cercetare (figura 8.1):

în Carpații Meridionali, în aproprierea localității Mălaia, județul Vâlcea;

în Carpații Meridionali, în apropierea Lacului de acumulare Vidra, județul Vâlcea;

Ca analiză preliminară a datelor, s-au determinat o serie de caracteristici de interes, prezentate în tabelul 9.1.

Tabelul 9.1. Caracteristici de interes pentru cele două seturi de date LiDAR

Analiza preliminară evidențiază diferențe semnificative între cele două suprafețe de studiu. Densitatea norului de puncte, unul dintre factorii care influențează acuratețea de reprezentare a suprafeței terestre, este semnificativ mai mare în cazul arealului Mălaia. Această diferență se explică prin faptul că acest areal a fost scanat în două campanii de zbor succesive (datele rezultate fiind suprapuse), în timp ce pentru zona Vidra s-a realizat o singură campanie de înregistrări.

Chiar dacă analiza preliminară ar justifica folosirea exclusivă a setului de date pentru arealul Mălaia, ca fiind superior calitativ, considerăm că este oportună cercetarea acurateții de modelare a suprafeței și pentru cazul mai puțin optim al arealului Vidra. În fapt, pentru zona de studiu Mălaia densitatea punctelor este ridicată pentru standardele actuale, astfel că analiza poate oferi o estimare a potențialului actual, chiar viitor, al tehnologiei.

9.3. Realizarea măsurătorilor geo-topografice

Scopul cercetării îl reprezintă estimarea potențialului de utilizare a tehnologiei LiDAR în silvicultură, ca sursă de date altimetrice pentru modelarea suprafeței terestre. Astfel, încă de la început s-a avut în vedere necesitatea obținerii unui set de date independent, folosit ca bază la estimarea gradului de acuratețe al datelor LiDAR. Un set de măsurători de control ar trebui să aibă un nivel de precizie și o densitate de eșantionare relativ apropiate față de tehnologia LiDAR. Dacă se iau în considerare condițiile de teren specifice scopului urmărit (acoperirea terenului cu vegetație forestieră), rezultă că singura sursă viabilă de date de control o reprezintă măsurătorile geo-topografice realizate cu aparatură modernă, de tip stație totală, și geo-referențiate folosind echipament de tip GNSS.

De la bun început, posibilitatea amplasării suprafețelor test în arealul Vidra a fost exclusă, datorită unui cumul de factori: relieful accidentat, densitatea foarte ridicată a arboretului care îngreunează utilizarea echipamentelor GNSS, precum și lipsa unei infrastructuri rutiere care să permită accesul cu aparatura geo-topografică, în condiții de siguranță.

Astfel, singura variantă fezabilă a reprezentat-o amplasarea suprafețelor de probă în cadrul arealului Mălaia. În etapa de planificare s-au stabilit, ca amplasare generală, un număr de opt suprafețe (figura 9.4), la delimitarea cărora s-a avut în vedere:

amplasarea în condiții de înclinare a terenului și consistență a vegetației cât mai variate;

situația infrastructurii rutiere, pentru a facilita accesul în condiții de siguranță dar și pentru a asigura aproprierea de zone deschise, necesare înregistrărilor GPS;

Măsurătorile topografice s-au desfășurat în condiții de lucru dificile, cu teren accidentat și vizibilitate îngreunată datorită frunzișului și trunchiurilor arborilor. Luând în seamă și numărul ridicat de puncte de radiat, randamentul redus de lucru a impus limitarea zonelor test la suprafețe de aproximativ 1-2000 metri pătrați.

Dintre cele opt suprafețe de probă, șapte au fost amplasate în interiorul pădurii (figura 9.1), cea de-a opta (identificată în continuare cu denumirea suprafață de referință) fiind amplasată în exteriorul acesteia, pentru control (figura 9.2).

Fig. 9.1. Suprafața de probă 1, amplasată în interiorul pădurii, în cadrul arealului Mălaia.

Cartarea suprafețelor de probă a presupus amplasarea a 27 puncte de stație, parte a opt drumuiri (cinci sprijinite și trei închise pe punctul de plecare). Pentru geo-referențierea măsurătorilor, coordonatele a 12 dintre punctele de stație (amplasate în zone deschise, în afara pădurii) au fost determinate cu echipamentul GPS Leica SR-20, folosind procedeul de măsurare „static cu post-procesare”. Această metodă presune folosirea simultană a cel puțin două receptoare GPS: unul dintre acestea (denumit bază) este instalat pe un punct cunoscut, înregistrând date pentru un anumit interval de timp. Un al doilea receptor GPS (denumit rover) înregistrează date pentru o serie de puncte noi. Atâta timp cât intervalul de înregistrare al rover-ului se încadrează în intervalul de timp în care receptorul bază este activ, rezultă prin post-procesare o serie de corecții aplicate măsurătorilor care sporesc precizia de poziționare.

În cazul de față, dintre cele 12 puncte poziționate prin GPS, 4 au fost folosite ca bază. S-a folosit o valoare pentru unghiul-mască de 15 grade (ignorarea semnalelor provenite de la sateliți aflați mai jos de 15 grade față de orizont) și un prag minim de 5 vectori (atâta timp cât nu sunt disponibile semnale de la minim 5 sateliți, înregistrarea datelor este întreruptă). Timpul de staționare pentru receptorul bază a fost setat la 90 minute (5400 epoci de măsurare), iar pentru echipamentul rover s-au staționat câte 30-45 minute (1800-2700 epoci) pe fiecare punct. Toate punctele determinate cu receptorul rover s-au aflat la cel mult un kilometru de receptorul bază la care au fost legate prin post-procesare. Astfel, majoritatea valorilor pentru indicatorul de precizie PDOP (Position Dilution of Precision) se situează în final între 1 și 3.

Fig. 9.2. Suprafața de probă pentru referință, amplasată în afara pădurii (parțial acoperită de o livadă), în cadrul arealului Mălaia.

Caracteristicile generale pentru cele opt suprafețe-test descrise anterior sunt prezentate în tabelul 9.2. Reprezentarea reliefului pentru suprafețele-test este inclusă în Anexa 1.

Tabelul 9.2. Caracteristici generale pentru suprafețele test delimitate în cadrul măsurătorilor topografice

Se observă faptul că, pentru o parte dintre suprafețele de probă, densitatea punctelor rezultate prin măsurători topografice este redusă prin comparație cu densitatea punctelor LiDAR din clasa Teren (valori de 0.1-0.2 față de 0.6-0.7). Trebuie amintit însă faptul că, în timp ce punctele LiDAR sunt distribuite aleatoriu, măsurătorile topografice sunt realizate manual, astfel că pot surprinde variația topografică prin mai puține puncte măsurate.

9.4. Filtrarea datelor LiDAR. Determinarea erorilor de filtrare.

9.4.1. Precizări generale

Filtrarea datelor LiDAR este procesul prin care ,din norul de puncte inițial, se separă acele puncte considerate a se afla la nivelul terenului. Aceste puncte, care sunt atribuite clasei denumită Teren, reprezintă ulterior datele de intrare pentru modelarea suprafeței terestre.

După cum s-a discutat anterior, filtrarea este prima etapă de prelucrare a datelor LiDAR asupra căreia s-a putut interveni în cercetarea de față. Astfel, unul dintre obiective l-a reprezentat determinarea efectului erorilor de filtrare asupra Modelului Digital Altimetric al Terenului (MDA).

9.4.2. Zona de studiu pentru analiza procesului de filtrare

Determinarea erorilor de filtrare s-a făcut în două variante:

prin comparare cu setul de LiDAR filtrat corect, referința fiind clasificarea pusă la dispoziție de către INCDS „Marin Drăcea”;

prin compararea cu o sursă de date independentă ( măsurători geo-topografice);

Pentru prima variantă s-a delimitat o suprafață extinsă, acoperind aproximativ 1.2 km2 (figura 9.3). Prin delimitarea suprafeței s-a urmărit surprinderea unui relief complex, cu expoziții variate și factor de fragmentare ridicat (figura 9.4). În cadrul acesteia, densitatea punctelor LiDAR este de 14.7 puncte/m2, densitatea punctelor din clasa Teren este de 1.7 puncte/m2, panta medie are valoarea de 32.60 grade iar gradul de acoperire cu vegetație este de 83 de procente. Pentru claritate, această suprafață va fi identificată în continuare ca suprafața de probă A.

Estimarea impactului filtrării asupra acurateții de reprezentare a suprafeței terestre s-a făcut deci pentru:

suprafața de probă A;

suprafețele de probă 1-8, măsurate în teren;

Menționăm că dintre cele opt suprafețe de probă pentru care s-au realizat măsurători topografice, suprafața nr. 3 este amplasată în interiorul suprafeței de probă A, iar restul se află în imediata vecinătate a acesteia (figura 9.3).

Fig. 9.3. Hartă cu amplasamentul pentru suprafețele pentru care s-au analizat erorile de filtrare a datelor LiDAR. Sistem de proiecție: UTM zona 35N.

Z„ Fig. 9.4. Model Digital Altimetric pentru suprafața de studiu A, aleasă pentru

analiza erorilor de filtrare (factor de exagerare verticală: 1.5).

Definirea erorii de filtrare

Filtrarea este procesul prin care se prelucrează un set de măsurători LiDAR, cu scopul de a separa din norul de puncte doar punctele aflate la nivelul terenului, ca etapă preliminară modelării suprafeței terestre. Filtrarea poate produce două tipuri de erori:

Erori de comitere: date de punctele care nu sunt eliminate, deși reprezintă obstacole aflate deasupra suprafeței terestre.

Erori de omitere: date de punctele aflate la nivelul terenului, eliminate în mod eronat de algoritmul de filtrare.

În fapt, nu interesează numărul punctelor incorect clasificate, ci efectul acestora asupra produsului final, adică modelul suprafeței terestre. Astfel, putem considera eroarea de filtrare ca o diferență de altitudine (eroare verticală) între suprafața obținută pe baza punctelor rezultate prin filtrare și o suprafață de referință. Pentru suprafața de probă A s-a folosit ca referință o suprafață interpolată pe baza clasificării manuale a punctelor din setul de date LiDAR, așa cum a fost pusă la dispoziție de către INCDS „Marin Drăcea”. Pentru suprafețele de probă nr. 1-8, suprafața de referință a fost obținută folosind punctele rezultate prin măsurători topografice.

Pentru determinarea erorilor verticale s-au parcurs următoarele etape:

filtrarea norului de puncte;

generarea prin interpolare a unui model raster pentru suprafața terenului (cu rezoluția de 0.5 metri) folosind ca date de intrare norul de puncte obținut în urma filtrării;

suprapunerea modelului raster cu un model al suprafeței de referință, obținut tot prin interpolare, pe baza datelor de control;

determinarea erorii verticale, pentru fiecare celulă a modelului, ca diferența dintre valorile de altitudine asociate celulei în cadrul celor două modele.

Rezultă astfel un model raster al erorilor, ale cărui celule se suprapun cu celulele suprafețelor interpolate. Pe baza acestui model, s-au calculat o serie de indicatori generali ai acurateții suprafeței:

eroarea medie, al cărei semn (pozitiv sau negativ) indică tendința de supraestimare sau subestimare a altitudinii;

eroarea medie absolută, indicator general al gradului de apropriere dintre suprafața obținută și suprafața de referință;

eroarea medie pătratică (RMSE – Root Mean Square Error), indicator global al calității produsului, din punct de vedere al acurateții.

Eroarea medie pătratică s-a calculat cu formula (Li, 1988):

[9.1]

unde

– RMSE este eroarea medie pătratică;

– E(i) erorile verticale asociate;

– n numărul de celule ale modelului raster;

Alegerea algoritmilor de filtrare

Până în prezent s-au dezvoltat numeroși algoritmi de filtrare, adaptați pentru diferite domenii de cercetare sau condiții de teren. Parcurgerea literaturii de specialitate reflectă faptul că nu există un consens în legătură cu gradul de acuratețe al algoritmilor, în general sau pentru situația terenurilor forestiere.

Astfel, s-a stabilit ca obiectiv secundar necesitatea alegerii unui grup de algoritmi care să fie supuși unei comparații riguroase, cantitative a performanței. La alegerea algoritmilor s-a urmărit ca aceștia să se încadreze în categorii diferite, ca metodă de filtrare. De asemenea, s-au ales atât algoritmi de filtrare generali, cât și algoritmi specifici pentru filtrarea terenului acoperit de vegetație forestieră. Un alt considerent a fost includerea, pe cât posibil, a algoritmilor disponibili în mod gratuit, în detrimentul programelor comerciale.

În final, a rezultat un set de nouă algoritmi de filtrare, șapte dintre care sunt disponibili gratuit. În funcție de modul de lucru, aceștia se încadrează în trei categorii principale (Cățeanu și Arcadie, 2017):

Algoritmi bazați pe interpolarea de suprafețe:

GroundFilter:

Algoritm bazat pe predicție liniară, implementat în pachetul software FUSION (Kraus și Pfeifer, 1998). Presupune interpolarea preliminară a unei suprafețe folosind întreg setul de date LiDAR, care se află deci între nivelul terenului și coronament. În continuare se atribuie ponderi fiecărui punct, în funcție de valoarea sa reziduală (distanța verticală până la suprafață). Principalii parametri prin care utilizatorul poate interveni asupra procesului sunt: valoarea de deplasare g (punctele cu un reziduu negativ mai mare decât g vor avea ponderea de 1.0) și parametrul w (punctele cu un reziduu pozitiv mai mare decât g+w vor avea o pondere de 0.0). Într-un proces iterativ, se vor genera suprafețe succesive, care sunt „atrase” spre punctele cu o pondere ridicată (care în principiu sunt mai aproape de suprafața reală a terenului). Numărul iterațiilor este ales de utilizator. Referirea la algoritm se va face în continuare prin termenul Groundfilter.

Height Filtering:

Algoritmul face parte din suita de programe BCAL (Boise Center Aerospace Laboratory) LiDAR Tools, fiind dezvoltat pentru filtrarea punctelor aflate la sol în cazul terenurilor acoperite cu vegetație de tip „pelin” (Streutker și Glenn, 2006). Se generează o rețea de celule, iar pentru fiecare celulă se identifică punctul de altitudine minimă. Pe baza acestui subset de puncte se generează o suprafață, care aproximează la nivel grosier suprafața terenului. Pentru toate punctele aflate sub suprafață, precum și cele aflate deasupra, dar până la o anumită distanță-limită, se atribuie clasa Teren. În continuare, algoritmul va fi identificat ca BCAL.

Multiscale Curvature Classification:

Distribuit gratuit sub formă de executabil, cu denumirea comercială MCC-LIDAR. Algoritmul realizează interpolarea unei suprafețe, la rezoluția definită de un parametru de scară. Se calculează curbura suprafeței, iar fiecare observație pentru care se depășește toleranța pentru curbură (stabilită de utilizator) este eliminată din setul de puncte. Aceste etape se repetă până se atinge un prag de convergență definit (procentul de puncte eliminate per iterație scade sub o anumită valoare). Procedura se urmează de trei ori, la domenii de scară multiple (0.5 x factor de scară, 1 x factor de scară, 1.5 x factor de scară). Toate punctele rămase în setul de date după convergența celui de-al treilea model sunt clasificate ca Teren. Algoritmul este dezvoltat pentru filtrarea datelor LiDAR în terenuri împădurite cu cantitate ridicată de biomasă (Evans și Hudak, 2007). În continuare, va fi identificat ca MCC.

gLidar:

Algoritmul este distribuit ca executabil gratuit de către GeMMA Lab. Folosind o abordare de sus-în-jos, o suprafață este interpolată și apoi condusă iterativ în jos, spre suprafața terenului. Pentru filtrare se calculează reziduurile pentru fiecare punct (distanța verticală dintre punct și suprafața curentă). Folosind o transformare matematică de tip top-hat, se compară reziduurile pentru grupuri de puncte învecinate, valorile extreme fiind identificate și eliminate (Mongus și Zalik, 2012). Pornind de la o rezoluție mai mare, se elimină succesiv din puncte, pentru rezoluții tot mai mici. În continuare, algoritmul va fi identificat ca gLiDAR.

Algoritmi care folosesc metoda TIN-adaptation (adaptarea Rețelei Neregulate de Triunghiuri)

Lasground:

Algoritmul este unul comercial, parte din suita de unelte pentru procesarea datelor LiDAR LASTools (Rapidlasso). Ca procedeu de filtrare, este similar cu algoritmul „TIN adaptiv” dezvoltat de către Axelsson (2000). Nefiind gratuit, informații detaliate despre modul de lucru sunt inaccesibile, dar la modul general un algoritm „TIN adaptiv” funcționează astfel: o Rețea Neregulată de Triunghiuri (TIN) este generată folosind un subset de puncte de pornire (seed points) de altitudine minimă. Numărul de triunghiuri al rețelei este sporit succesiv, prin adăugare de puncte din setul inițial de date. La fiecare iterație, se verifică distanța și unghiul fiecărui punct până la fațeta rețelei aflată sub acesta. Când valorile sunt sub un prag-limită, punctul este adăugat la rețea. Procesul se încheie când niciunul din punctele aflate în afara rețelei nu îndeplinește condițiile de adăugare la rețea. Modelul TIN final aproximează suprafața terenului, iar nodurile sale (care sunt puncte din setul de date inițial) sunt atribuite clasei Teren. Implementarea din programul LASTools folosește un parametru denumit „extra-fine” care intensifică procesul de identificare a punctelor de pornire, pentru a îmbunătăți rezultatul în cazul terenurilor accidentate.

Lasground-new:

Algoritmul Lasground-new este o versiune mai nouă a algoritmul Lasground descris anterior. Scopul algoritmului este îmbunătățirea procesului de filtrare pentru terenurile complexe (de exemplu dealuri abrupte învecinate cu zone construite). Este de asemenea inclus în pachetul de unelte LASTools.

Algoritmi bazați pe morfologia matematică

Elevation Threshold with Expand Filter:

Algoritm dezvoltat de Zhang și Whitman (2005), de asemenea implementat în pachetul ALDPAT. Presupune generarea unei rețele de celule care este suprapusă cu datele LiDAR, punctul cu cea mai mică altitudine fiind identificat pentru fiecare celulă. Toate celelalte puncte sunt considerate a se aflate deasupra nivelului terenului și deci eliminate. În iterația următoare, dimensiunea celulelor este dublată. În interiorul fiecărei celule, pentru fiecare punct se calculează diferența de nivel față de punctul de altitudine minimă. Acele puncte pentru care diferența de nivel depășește o valoare-limită sunt eliminate. Pragul este determinat ca produsul dintre un parametru de înclinare (setat de utilizator) și dimensiune curentă a celulei. Procesul se încheie atunci când printr-o iterație nu se elimină niciun punct, sau când se atinge numărul de iterații stabilit de utilizator. În continuare, algoritmul va fi identificat ca ETEW.

Maximum Local Slope:

Algoritm dezvoltat de către Zhang și Whitman (2005) și implementat în pachetul gratuit ALDPAT (Airborne LiDAR Data Processing and Analysis Tools). Algoritmul generează o rețea de celule (cu dimensiunea celului setată de utilizator) pe care o suprapune pe setul de date LiDAR. Pentru fiecare punct se determină înclinarea față de toate punctele învecinate (aflate în interiorul unei zone circulare de rază aleasă de utilizator). Dacă oricare dintre valorile de înclinare depășește un prag stabilit de utilizator (parametrul slope), observația va fi încadrată în clasa Teren. În continuare, algoritmul va fi identificat ca MLS.

Simple Morphological Filter:

Algoritmul este disponibil gratuit ca unealtă pentru mediul de programare MATLAB (MathWorks©). Ca principiu, se folosește morfologia matematică pentru identificarea punctelor aflate la nivelul terenului. Inițial se generează o suprafață de altitudine minimă, folosind cel mai jos punct din fiecare celulă a unei rețele rectangulare. Apoi, o secvență iterativă de operațiuni morfologice (numite „deschideri” și „închideri”) sunt realizate, ținând cont de doi parametri setați de utilizator: mărimea ferestrei (depinde de dimensiunea maximă a obiectelor aflate deasupra terenului) și înclinarea maximă a terenului. În urma acestor operațiuni, fiecare celulă va fi etichetată ca Obstacol sau Teren. Elevația terenului corespondentă celulelor Obstacol se estimează prin interpolarea valorilor pentru celulele Teren învecinate. Astfel, suprafața reprezintă un model grosier al terenului, folosit la clasificarea punctelor. Atunci când diferența de nivel dintre un punct și model se află sub o valoare-prag (calculată ca suma dintre o toleranță constantă și un factor de scalare aplicat înclinării terenului), acel punct este atribuit clasei Teren (Pingel et al, 2013). În continuare, algoritmul va fi identificat ca SMRF.

Optimizarea procesului de filtrare

Fiecare dintre algoritmii de filtrare aleși pentru analiza acurateții implică setarea unor parametri de către utilizator, prin care se controlează procesul de filtrare. Acolo unde prin studierea literaturii de specialitate s-au găsit menționate valorile optime pentru parametri, acestea au folosite. Unde acest lucru nu a fost posibil, s-a apelat la următoarea procedură:

pentru fiecare algoritm, se pornește de la valorile standard ale parametrilor;

valoarea fiecărui parametru se modifică (în plus sau minus), rezultatele fiind analizate vizual;

dacă modificarea parametrului nu afectează rezultatul filtrării în mod vizibil, parametrul se elimină din analiză, rămânând la valoarea standard;

dacă modificarea parametrului are efect vizibil asupra filtrării, se continuă modificarea acestuia până când valoarea nu mai apar modificări vizibile – rezultă astfel o listă de valori de testat pentru fiecare parametru;

pentru a testa interacțiunea dintre parametri, se realizează combinațiile dintre valorile alese ca fiind reprezentative pentru fiecare parametru.

Prin procedura descrisă mai sus, a rezultat un număr de combinații ale valorilor parametrice pentru fiecare algoritm, situat între 25 (pentru algoritmul ETEW) și 104 (pentru algoritmul Groundfilter). În total, pentru cei nouă algoritmi s-au analizat 568 de combinații.

Analizarea unei combinații de valori parametrice presupune următoarele etape:

filtrarea datelor LiDAR pentru suprafața de probă A, folosind setul de valori parametrice pentru combinația curentă;

generarea unui model al suprafeței cu rezoluția de 1.0 metri, prin interpolare cu algoritmul IDW (Inverse Distance Weighted), folosind setul de puncte filtrate;

calcularea erorii medie pătratice a modelului (RMSE), pe baza erorilor verticale ale celulelor (diferența de altitudine dintre modelul interpolat și suprafața de referință);

pentru fiecare algoritm, se identifică acea combinație de parametri care produce cea mai precisă modelare a suprafeței terestre (cea mai mică eroare medie pătratică).

Rezultă astfel nouă seturi de puncte filtrate, câte unul pentru fiecare dintre algoritmii testați. Aceste seturi de date sunt folosite în continuare pentru analiza comparativă între algoritmi. Valorile parametrilor identificate ca optime urmând procedura de mai sus sunt prezentate în tabelul 9.3.

Tabelul 9.3. Valori finale pentru parametrii algoritmilor de filtrare a datelor LiDAR; doar acei parametri care au alte valori decât cele standard sunt prezentați.

Distribuția spațială a erorilor de filtrare

Norul de puncte LiDAR, ca obiect asupra căruia se aplică procedura de filtrare a observațiilor, este orientat în spațiu – fiecărui punct i se asociază un set de coordonate X, Y, Z. Astfel, filtrarea datelor LiDAR este un proces cu o componentă spațială. Este de interes deci analizarea erorilor de filtrare nu doar ca distribuție spațială, ci și ca distribuție în spațiu. Pentru a reduce din variabilitatea erorilor de filtrare, cu scopul de a simplifica analiza, s-au delimitat o serie de clase ale erorilor verticale, astfel:

erori foarte reduse, între 0.0 și 0.2 metri;

erori reduse, între 0.2 și 0.5 metri;

erori semnificative, între 0.5 și 1.0 metri;

erori majore, peste 1.0 metri.

Facem mențiunea că delimitarea este una subiectivă, întrucât nu există un standard în domeniu privind acest aspect. S-a ținut totuși seama de gradul de acuratețe urmărit în mod uzual în domeniul forestier, precum și de faptul că în practică modelarea terenului s-ar face probabil la nivel de parchet. Din aceste considerente, considerăm spre exemplu o eroare de determinare a altitudinii de până la 0.2 metri ca fiind foarte redusă, prin prisma efectului acesteia asupra soluțiilor practice adoptate.

După clasificarea erorilor, s-au generat hărțile cu distribuția spațială a claselor, calculându-se de asemenea suprafața acoperită de fiecare dintre acestea.

De asemenea, s-a analizat gradul de auto-corelație spațială a erorilor de filtrare, sau măsura în care erorile tind să se organizeze în grupuri de valori apropiate. Cu alte cuvinte, auto-corelația spațială determină gradul în care o observație influențează sau depinde de observațiile aflate în apropriere. O grupare a unui număr de valori apropiate, care depășește probabilitatea de grupare aleatorie poartă numele de cluster (Webster și Oliver, 2007). Gruparea unor valori ridicate (relativ la media setului de observații) se numește cluster pozitiv, iar gruparea valorilor reduse (relativ la medie) se numește cluster negativ.

Pentru cuantificarea gradului de auto-corelație spațială, s-a folosit tehnica LISA propusă de Anselin (1995). Unul dintre indicatorii propuși de autor pentru estimarea auto-corelației spațiale este Local Moran’s I, calculat cu formula:

[9.2]

unde:

I este indicatorul Moran’s I;

Zi deviația observației i, relativ la medie;

Wij matrice de ponderi ;

N numărul total de observații;

S0 suma algebrică a valorilor matricei Wij;

Indicele Local Moran’s I are valori între -1 și +1, unde:

-1 reprezintă o organizare perfectă a valorilor nesimilare (dispersie perfectă);

0 lipsa corelației (distribuție perfect aleatorie)

+1 organizare perfectă pe clustere (opusul dispersiei).

Analiza distribuției spațiale a erorilor s-a efectuat pentru suprafața de probă A, în cazul filtrării cu ajutorul algoritmilor: Lasground-new, MLS și gLiDAR.

Factori de influență a erorilor de filtrare

Este de așteptat ca erorile de filtrare să fie influențate de condițiile de teren și de caracteristicile datelor LiDAR. Factorii care au fost identificați ca având o potențială legătură cu amplitudinea erorilor sunt:

Înclinarea terenului, care s-a determinat pe baza Modelului Digital Altimetric de referință.

Rugozitatea suprafeței terenului, calculată pe baza Indicelui de Rugozitate al Terenului (TRI).

Indicele de Rugozitate al Terenului cuantifică gradul de omogenitate topografică a suprafeței terestre, prin prisma diferențelor de nivel dintre celulele Modelului Digital Altimetric (Riley ș.a., 1999). Calcularea indicelui se face pentru fiecare celulă a modelului în parte, cu formula:

[9.3]

unde:

TRI este Indicele de Rugozitate al Terenului;

Z0 altitudinea celulei curente;

Zij altitudinea pentru cele 8 celule învecinate celulei curente;

Curbura maximă a suprafeței, este valoarea maximă dintre curbura în profil și curbura în plan a suprafeței.

Densitatea punctelor LiDAR, care s-a determinat sub forma unei rețele de celule rectangulare, cu suprafața de 1 m2, pentru fiecare celulă fiind calculat numărul de puncte aflate în interiorul acesteia. Pentru reducerea nivelului de zgomot al modelului, s-a aplicat un filtru de medie cu dimensiunea ferestrei de 5 x 5 celule.

Gradul de acoperire cu vegetație, pentru estimarea căreia s-au folosit acele observații din setul de date LiDAR care sunt atribuite altor clase decât Teren (de regulă Vegetație joasă, Vegetație medie sau Vegetație înaltă). Determinarea densității s-a făcut într-o rețea rectangulară de celule, cu formula (McGaughey, 2004):

[9.4]

unde:

– δi este densitatea coronamentului (exprimată în procente) pentru celula i;

– n numărul de puncte din celula i situate deasupra unui prag de înălțime ales (în cazul de față, doi metri deasupra terenului);

– N numărul total de puncte din interiorul celulei i;

Limitele claselor sunt prezentate în tabelul 9.4. În lipsa unor valori standardizate, clasele s-au delimitat în funcție de distribuția statistică a valorilor. Semnificația atribuită claselor este una descriptivă, relativă la suprafața aleasă pentru analiza procesului de filtrare (suprafața de probă Mălaia).

Distribuția spațială a factorilor analizați, reprezentată grafic sub formă de hărți tematice, se regăsește în Anexa 2.

Pentru analizarea efectului acestor factori asupra erorilor, s-au ales trei dintre algoritmii de filtrare ca fiind reprezentativi, câte unul din fiecare categorie. Aceștia sunt:

gLiDAR, pentru clasa algoritmilor bazați pe interpolarea de suprafețe;

Lasground-new, pentru clasa de algoritmi bazați pe adaptarea Rețelei Neregulate de Triunghiuri (TIN);

MLS, pentru algoritmii de tip morfologic;

Tabelul 9.4. Gruparea pe clase a factorilor luați în considerare la analiza erorilor de filtrare.

Interpolarea suprafeței terenului. Determinarea erorilor de interpolare.

Precizări generale

Interpolarea este procesul prin care un set de observații punctuale este folosit la generarea unei suprafețe continue, care exprimă variabilitatea atributului studiat. În cazul de față, observațiile sunt reprezentate de punctele obținute prin tehnologia de scanare laser, pentru care se cunoaște poziția planimetrică (coordonatele X și Y) și altimetrică (coordonata Z). Folosind ca date de intrare punctele din clasa Teren, prin interpolarea valorilor Z va rezulta o suprafață continuă sub forma unui structuri raster (celule rectangulare organizate pe rânduri și coloane). Această reprezentare este Modelul Digital Altimetric (MDA).

Interpolarea este una din etapele din procesul de procesare a datelor LiDAR care va produce erori, influențând acuratețea de reprezentare a suprafeței terestre. Astfel, unul dintre obiectivele cercetării a fost determinarea efectului erorilor de interpolare asupra produsului final (Modelul Digital Altimetric).

Zona de studiu pentru analiza procesului de interpolare

Pentru estimarea și analizarea erorilor de filtrare, s-a folosit setul de date LiDAR pentru arealul Vidra. Întrucât acesta are o formă neregulată (figura 8.3) care poate influența negativ procesul de interpolare, s-a delimitat o zonă sub formă regulată de paralelogram, în interiorul acestui areal (figura 9.5). Totodată, prin această delimitare s-a eliminat secțiunea sudică din setul de date, lipsită de vegetație forestieră. A rezultat suprafață de probă Vidra, care acoperă 160 de hectare, pe care vegetația forestieră este prezentă aproape în întregime. Altitudinea variază între 1350 și 1750 metri și înclinare medie a terenului este de 26 de grade. Densitatea medie a observațiilor LiDAR, reprezentată de numărul de puncte pe unitatea de suprafață, este de 5.2 puncte/m2, respectiv 0.89 puncte/m2 pentru punctele atribuite clasei Teren.

Date pentru validarea interpolării

După cum s-a menționat anterior, din considerente de timp, costuri sau siguranță, nu s-au putut efectua măsurători topografice în arealul Vidra. Astfel, fără o sursă de date independente, care să poată servi ca referință pentru estimarea erorilor de interpolare, s-a recurs la metoda numită cros-validare (Davis, 1987). Tehnica este folosită în mod uzual în studiile de analiză a erorilor de măsurare, în situația în care nu sunt disponibile măsurători independente pentru validare (Anderson et al, 2005; Aguilar, 2005).

Fig. 9.5. Suprafața de probă Vidra , folosită pentru analiza procesului de interpolare.

În esență, metoda presupune împărțirea, în mod aleatoriu, a setului de observații în două subseturi: unul pentru predicție (care cuprinde majoritatea observațiilor) și unul pentru validare (care cuprinde o mică parte dintre observații). Pentru cercetarea de față, s-a optat pentru o împărțire 95 procente observații pentru predicție cu 5 procente observații pentru validare. Proporția a fost aleasă pentru ca testarea acurateții să fie semnificativă, dar fără a se compromite integritatea datelor de predicție (Bater & Coops, 2009).

Astfel, din setul de observații LiDAR inițial, s-a extras aleatoriu un număr reprezentând 5 procente din total (n = 52358), care au servit la validare. Restul de 95 de procente din observații (aprox. 1.06 milioane puncte) au fost folosite pentru realizarea interpolărilor.

Definirea erorii de interpolare

Interpolarea urmărește generarea unei suprafețe continue pornind de la un set de observații discrete. Astfel, este necesară estimarea valorilor pentru locații unde nu există observații, deci procesul este afectat de un anumit grad de incertitudine. Această incertitudine se poate cuantifica prin determinarea și analizarea erorilor de interpolare.

Pentru aceasta, suprafața interpolată se suprapune cu observațiile de validare. Pentru fiecare din punctele de validare se cunosc două valori: cea specifică punctului (considerată corectă) și cea interpolată (considerată incertă). Astfel, pentru fiecare punct de validare se poate calcula eroarea de interpolare, cu formula:

[9.5]

unde:

E(i) este eroarea verticală pentru punctul i;

ZDTM(i) valoarea interpolată pentru poziția planimetrică a punctului i;

ZLiDAR(i) valoarea de altitudine a punctului i, așa cum a rezultat în urma măsurătorilor;

Odată calculate erorile verticale individuale pentru cele n observații din setul de validare, o serie de indicatori statistici oferă informații privind acuratețea generală a suprafeței interpolate:

eroarea medie, al cărei semn (pozitiv sau negativ) indică prezența unei tendințe de supraestimare sau subestimare a valorilor de elevație;

eroarea medie absolută, indicator general al gradului de apropriere dintre suprafața obținută și suprafața de referință;

eroarea medie pătratică (RMSE – Root Mean Square Error), calculată cu formula 9.1.

Alegerea algoritmilor de interpolare

Până în prezent, s-au dezvoltat numeroase metode de interpolare, menite să estimeze valori ale atributelor analizate pentru locații ne-eșantionate (Mitas și Mitasova, 1999). Astfel, s-a remarcat necesitatea de a testa mai mulți algoritmi de interpolare, pentru a realiza analiza comparativă a gradului de acuratețe oferit de aceștia. La alegerea algoritmilor s-a urmărit ca aceștia să fie de diferite tipuri (conform clasificării descrise în secțiunea 6.2.). De asemenea, se are în vedere că același algoritm, implementat în programe software diferite, nu produce în mod necesar rezultate identice (Pietrzak, 2013, citat în Sterenczak, 2016). Din acest motiv, s-a urmărit ca algoritmii selectați să fie disponibili într-un pachet software comun.

În final s-a selectat un număr de nouă interpolatori, integrați în pachetul software de analiză și modelare spațială SAGA GIS. O parte dintre algoritmi necesită setarea unor valori parametrice de către utilizator; valorile testate sunt prezentate în tabelul 9.5. Procedura de interpolare specifică celor nouă algoritmi este detaliată în continuare:

Inverse Distance Weighted (IDW)

IDW este un interpolator deterministic local, de tip convex. Valorile exprimă rezultă folosind ecuația propusă de Shepard (1968):

[9.6]

unde:

Z este valoarea prezisă pentru punctul necunoscut x0 ;

w(di) funcția de ponderare a celor n puncte învecinate punctului x0 ;

z(xi) valoarea măsurată pentru punctul xi ;

di distanța dintre punctele x0 și xi ;

Funcția de ponderare reflectă principiul auto-corelării spațiale (valorile apropiate tind să fie asemănătoare) și are forma generală:

[9.7]

unde:

w este ponderea atribuită punctului xi la predicția valorii pentru o locație necunoscută x0 ;

di distanța dintre punctele x0 și xi ;

Atunci când se dorește ca punctele îndepărtate (față de locația pentru care se realizează predicția) să aibă o influență redusă, se atribuie o valoare mai mare exponentului u (Weng, 2006). Valoarea standard pentru parametrul u este 2 (caz în care ponderea scade cu pătratul distanței), deși în anumite situații valori precum 1 sau 3 conduc la o predicție mai precisă (Aguilar, 2005).

În cazul de față s-au testat patru valori pentru exponentul u: 0.5, 1, 2, 3.

Nearest Neighbour (NeN)

Algoritmul Nearest Neighbour este un algoritm deterministic local de tip convex, care presupune generarea de poligoane Thiessen. Poligonul Thiessen asociat unui punct reprezintă aria de influență a acelui punct, și se construiește unind bisectoarele perpendiculare dintre vecinii punctului (Brassel și Reif, 1979). Fiecare punct cunoscut se află în centrul unui poligon, iar toate punctele din poligon sunt mai apropiate de acel punct cunoscut decât de oricare altul (Isaaks și Srivastava, 1989). Astfel, suprafeței închise de un poligon Thiessen i se va atribui valoarea punctului cunoscut din centrul poligonului.

Natural Neighbour (NN)

Algoritmul Natural Neighbour este un interpolator deterministic local, de tip convex, dezvoltat de Sibson (1981). La fel ca în cazul algoritmului Nearest Neighbor, se construiește rețeaua de poligoane Thiessen. Apoi se inserează în rețea punctele pentru care se realizează predicția, generându-se noi poligoane Thiessen. Aceste poligoane se vor suprapune cu cele inițiale; fiecare poligon nou preia o parte din suprafața unui număr de poligoane inițiale (Sibson denumește punctele cunoscute din centrul acestor poligoane “vecinii naturali” ai punctului nou). Poligonul Thiessen al punctului nou este compus din părți ale poligoanelor vecinilor săi naturali. Ponderea atribuită punctelor cunoscute va depinde de proporția în care poligonul Thiessen asociat acoperă suprafața poligonului pentru care se face predicția.

Delauney Triangulation (DT)

Algoritmul de interpolare Delauney Triangulation generează o Rețea Neregulată de Triunghiuri (TIN) pentru setul de observații cunoscute. Rezultă un model vectorial care unește punctele în spațiul tridimensional. Orice punct nou pentru care se estimează valoarea se va afla în interioul unui triunghi al rețelei. Valoarea prezisă pentru acel punct va depinde de valorile vârfurilor triunghiului respectiv, fiind calculată printr-o interpolare liniară sau polinomială (Ripley, 1981).

Ordinary Kriging (OK)

Algoritmul Ordinary Kriging este un interpolator geostatistic care urmărește cuantificarea gradului de autocorelație spațială a observațiilor. Aceasta se face printr-un grafic numit semi-variogramă, care rezumă forma de ansamblu a variației (figura 9.6), precum și magnitudinea și scala spațială a acesteia (Oliver & Webster, 1990). Pe acest grafic se potrivește un model de semivariogramă, care poate avea diferite forme: polinomial, circular, sferic, exponențial, gaussian etc. Funcția modelului este folosită pentru a atribui ponderi punctelor cunoscute. Ponderile depind de distribuția spațială a valorilor măsurate și sunt deci o funcție a co-varianței spațiale (Weng, 2006).

Tabelul 9.5. Caracteristicile algoritmilor de interpolare analizați și parametrii setați.

Fig. 9.6. Model de semi-variogramă pentru variabila Z.

Cu linie roșie se reprezintă funcția care modelează cel mai bine variația lui Z

(în acest caz, o funcție polinomială de gradul 2).

Radial Basis Functions (RBF)

Cu denumirea de Radial Basic Functions se identifică o clasă de algoritmi de interpolare exacți, care urmăresc potrivirea unei suprafețe flexibile (de tip Spline) printr-un set de puncte cunoscute (Erdogan, 2009). Funcția de interpolare trece prin punctele cunoscute (sau cât mai aproape de acestea), dar în același timp urmărește generarea unei suprafețe de tensiune minimă (Mitas & Mitasova, 1999). Algoritmii RBF pornesc de la presupunerea că există o anumită eroare de măsurare a punctelor cunoscute, care va fi redusă prin nivelarea la nivel local a suprafeței (Hengl, 2009). Pentru a realiza această nivelare locală se reduce la minim tensiunea suprafeței generate.

Există mai mulți algoritmi de tip RBF, care folosesc diferite funcții matematice pentru a defini suprafața. Patru astfel de algoritmi au fost incluși în cercetarea de față:

Algoritmul Multilevel B-Spline (BS) folosește o ierarhie de matrici de control pentru a genera o secvență de funcții Spline a căror sumă tinde la funcția de interpolare dorită (Lee ș.a., 1997).

S-au testat trei valori pentru numărul de matrici (sau nivele) generate de algoritm: 5, 10, 14.

Algoritmul Cubic Spline (CS) generează suprafața printr-o funcție continuă de tip Spline cubic bivariat (Haber ș.a., 2001). Toleranța setată de utilizator are efect asupra senzitivității funcției spline obținute. Valori mai mici ale parametrului conduc la o reprezentare a suprafeței mai puțin variabilă la nivel local.

S-au testat trei valori pentru toleranță: 80, 140, 200.

Algoritmul Thin Plate Spline (TPS) folosește funcția cu același nume pentru generarea suprafeței. Funcția Thin Plate Spline reprezintă generalizarea în spațiul bidimensional a funcției Spline Cubic.

Algoritmul Thin Plate Spline with TIN (TPSTIN) este o variantă a TPS care presupune generarea unei rețele TIN, anterior interpolării. În acest caz, în locul unei funcții TPS globale, se calculează câte o funcție TPS pentru fiecare triunghi al rețelei. Timpul de procesare este deci mult mai ridicat.

Pentru a reduce efectul de margine, la toate interpolările a fost folosită o zonă-tampon (buffer) la distanța de 5 metri.

Alegerea rezoluției spațiale pentru modelul suprafeței

Rezoluția modelului de suprafață (dată de dimensiunea celulelor rețelei) este unul dintre factori care influențează acuratețea reprezentării. Din acest motiv, s-a considerat că este oportună analizarea modelelor la diferite rezoluții spațiale.

Hengl (2006) propune o serie de reguli analitice pentru determinarea rezoluției optime de interpolare, pornind de la valori cunoscute precum: scara la care s-a realizat eșantionarea, dimensiunile și forma celor mai mici obiecte ce se doresc a fi reprezentate sau gradul de auto-corelare spațială.

Metodele descrise de autor fac referire la frecvența Nyquist propusă de Shannon (1949), care teoretizează că un semnal poate fi reconstruit dacă frecvența de sondare este de două ori mai mare decât frecvența de transmisie. Prin analogie, un model al suprafeței care să conserve cea mai mare parte din informația setului de date inițial ar avea o mărime a celulei egală cu jumătate din media spațierii minime dintre perechile de puncte (Bater & Coops, 2009).

În cazul de față, media distanțelor minime dintre puncte are valoarea de 0.56 metri. Întrucât testarea unei rezoluții spațiale de 0.25-0.30 metri implică resurse de procesare semnificative, s-au testat următoarele rezoluții: 0.5, 1.0, 1.5 și 2.0 metri.

Astfel pentru cei nouă algoritmi de interpolare, unii dintre ei cu mai multe variante de setare a parametrilor, au rezultat în final 92 de modele raster de analizat (câte 23 de modele pentru fiecare rezoluție).

Factori de influență a erorilor de interpolare

Interpolarea este un proces inexact, care presupune estimarea valorilor pentru locații noi, ne-eșantionate. Din acest motiv, procesul va fi afectat de un anumit grad de incertitudine, adică de erori de interpolare. Analiza literaturii de specialitate din domeniu a scos la iveală o serie de factori externi care pot avea efect asupra mărimii erorilor de interpolare. O parte dintre aceștia au putut fi calculați folosind datele avute la dispoziție, astfel că au fost incluși în analiză.

La analiza erorilor asociate procesului de interpolare s-au luat în considerare aceiași factori identificați pentru procesul de filtrare a datelor LiDAR:

Înclinarea terenului.

Indicele de Rugozitate al Suprafeței (TRI – Terrain Rugedness Index)

Curbura maximă a suprafeței.

Gradul de acoperire cu vegetație.

Densitatea punctelor de predicție.

Descrierea factorilor menționați mai sus se regăsește în secțiunea 9.4.7.

Pentru analiza stratificată a erorilor, factorii amintiți mai sus au fost grupați pe clase. Limitele claselor pentru fiecare factor sunt prezentate în tabelul 9.6. În lipsa unor valori standardizate, clasele s-au delimitat în funcție de distribuția statistică a valorilor. Semnificația atribuită claselor este una descriptivă, specifică suprafeței pentru care s-a analizat procesul de interpolare (arealul Vidra).

Distribuția spațială a factorilor analizați, reprezentată grafic sub formă de hărți tematice, se regăsește în Anexa 3.

Tabelul 9.6. Gruparea pe clase a factorilor luați în considerare la analiza erorilor de interpolare.

10. REZULTATE ȘI DISCUȚII

10.1. Filtrarea datelor LiDAR.

10.1.1. Precizări generale

Filtrarea datelor LiDAR, identificată ca una dintre sursele de erori care influențează calitatea Modelului Digital Altimetric, a fost analizată prin testarea a nouă algoritmi de filtrare (enumerați în tabelul 9.3). Procesul de filtrare s-a realizat în două variante (secțiunea 9.4.2):

Pentru o suprafață extinsă, de aprox. 1.2 km2, identificată ca suprafața de probă A din cadrul arealului Mălaia (figura 9.3); evaluarea acurateții de filtrare s-a făcut prin compararea cu datele LiDAR clasificate corect, așa cum au fost puse la dispoziție de INCDS „Marin Drăcea”.

Pentru un număr de opt suprafețe de probă, pe care s-au efectuat măsurători geo-topografice, identificate ca suprafețe de probă nr. 1-8 (figura 9.4); evaluarea acurateții de filtrare s-a făcut în acest caz prin compararea cu măsurătorile geo-topografice, ca sursă de date independentă.

Pentru fiecare situație analizată s-au calculat valorile pentru eroarea medie și pentru eroarea medie pătratică (RMSE) – acestea sunt prezentate în tabelul 10.1.

În primul rând, pentru suprafața de probă nr. 8 (care a fost delimitată în afara domeniului forestier, pentru a servi ca referință – figura 9.2) se remarcă faptul că valorile pentru RMSE și eroarea medie relativă sunt similare pentru majoritatea algoritmilor de filtrare. Pe această suprafață terenul este în mare parte descoperit, astfel că filtrarea norului de observații LiDAR se desfășoară în condiții optime. Astfel, șapte dintre algoritmi (MCC, gLiDAR, Lasground, Lasground-new, ETEW, MLS și SMRF) au valori ale RMSE între 0.25 și 0.27 metri, în timp ce valorile erorii medii sunt de 0.10-0.11 metri, indicând o supraestimare minoră a altitudinii. Ceilalți doi algoritmi de filtrare, Groundfilter și BCAL au o performanță mult mai scăzută, valorile indicelui RMSE asociate acestora fiind de 0.50 și respectiv 0.49 metri. Cei doi algoritmi au o ușoară tendință de subestimare a altitudinii, având valori ale erorii medii de -0.10 și respectiv -0.11 metri.

Se poate deci afirma că o valoare de 0.25-0.30 metri pentru eroarea medie pătratică reprezintă o referință față de care se pot raporta celelalte rezultate, pentru suprafețele în care condițiile de filtrare nu sunt optime (datorită prezenței vegetației forestiere).

În ceea ce privește suprafețele de probă nr. 1-7, amplasate în interiorul pădurii, se observă diferențe semnificative ale erorilor de filtrare. Eroarea medie pătratică variază de la o suprafață de probă la alta, chiar dacă algoritmul rămâne constant. Spre exemplu, algoritmul SMRF are valori ale erorii medii pătratice între 0.32 metri pentru suprafața de probă 6 și 1.30 metri pentru suprafața de probă 3. Diferențe apar și între algoritmi, la filtrarea punctelor pentru aceeași suprafață de probă. De exemplu, pentru suprafața de probă 6 eroarea medie pătratică are valori între 0.31 și 2.02 metri.

Pentru suprafața de probă A, indicele RMSE are valori între 0.34 și 2.25 metri, în timp ce eroarea medie relativă indică o supraestimare a altitudinii terenului pentru toți algoritmii în afară de BCAL. Analiza erorilor asociate acestei suprafețe confirmă observațiile pentru suprafețele de probă 1-7, care au o întindere mai redusă. Astfel, se remarcă o performanță relativ scăzută a algoritmilor BCAL și Groundfilter, în timp ce între ceilalți algoritmi performanța este în general apropiată.

Algoritmii de filtrare care au fost testați se împart în trei categorii principale: algoritmi bazați pe interpolare (Groundfilter, BCAL, MCC și gLiDAR), bazați pe adaptarea Rețelei Neregulate de Triunghiuri (Lasground și Lasground-new) și algoritmi morfologici (ETEW, MLS și SMRF). Calcularea erorilor pătratice medii între algoritmi, pe clase, arată faptul că algoritmii morfologici și cei care folosesc metoda adaptare-TIN produc erori de filtrare semnificativ reduse, prin comparație cu algoritmii care folosesc interpolarea suprafețelor pentru filtrare (figura 10.1). Acuratețea ridicată oferită de algoritmii de filtrare de tip adaptare-TIN este de asemenea remarcată de cercetări anterioare (Korzeniowska ș.a., 2014; Montealegre ș.a., 2015).

Fig. 10.1. Acuratețea medie de filtrare, pentru cele trei clase de algoritmi analizați.

Pentru fiecare categorie de algoritmi, acuratețea procesului de filtrare este mai ridicată pentru suprafața de probă A decât pentru suprafețele nr. 1-7. Posibile explicații sunt: gradul foarte ridicat de acoperire cu vegetație specific suprafețelor nr. 1-7 (care limitează penetrarea la sol a undelor laser), influența erorilor asociate etapei de măsurători terestre și efectul procesului de interpolare a suprafețelor. De asemenea, trebuie ținut seamă de faptul că algoritmii de filtrare au fost optimizați pentru suprafața de probă A, aceleași valori parametrice fiind apoi folosite și pentru suprafețele nr. 1-8 (vezi secțiunea 9.4.5).

O reprezentare vizuală a rezultatului filtrării se regăsește în Anexa 4.

10.1.2. Numărul de puncte rezultat în urma filtrării

Prin procesul de filtrare se identifică acele observații din norul de puncte LiDAR care este probabil să se afle la nivelul solului. Acestor puncte li se atribuie clasa Teren, restul punctelor fiind eliminate din setul de date. Pentru suprafața de probă A, numărul de puncte filtrare variază între 0.59 și 3.61 milioane de puncte în funcție de algoritm (tabelul 10.2). Numărul de puncte atribuite clasei Teren, conform clasificării puse la dispoziție de INCDS „Marin Drăcea”, este de 1.3 milioane.

Se constată că doar doi algoritmi filtrează un număr de puncte mai mic decât cel corect. Aceștia sunt:

Groundfilter, care filtrează 45 la sută din numărul corect de puncte și are o acuratețe scăzută relativ la ceilalți algoritmi (eroare medie pătratică de 1.23 metri pentru suprafața de probă A.

MLS, care filtrează 69 la sută din numărul corect de puncte și are o acuratețe bună raport la ceilalți algoritmi (eroare medie pătratică de 0.56 metri pentru suprafața de probă A).

Pentru ceilalți algoritmi, numărul de puncte îl depășește pe cel real, de până la aproape 3 ori (în cazul SMRF, cu 3.6 milioane puncte filtrate). Rezultă deci că majoritatea algoritmilor sunt afectați de erori de comitere mai semnificative decât erorile de omitere.

Cu alte cuvinte, este mai probabil ca un punct reprezentând un obstacol să fie etichetat incorect ca reprezentând suprafața terenului, decât invers. Întrucât valorile pentru RMSE sunt relativ bune, se presupune că o parte semnificativă din aceste puncte incorect atribuite clasei Teren se află la joasă înălțime față de sol. Astfel impactul lor asupra acurateții de modelare a suprafeței terestre este relativ redus. De exemplu, deși algoritmul Lasground-new filtrează un număr de puncte aproape dublu față de cel real (2.37 față de 1.3 milioane), valoarea RMSE asociată acestuia este de 0.34 metri (pentru suprafața de probă A).

În ceea ce privește suprafețele de probă nr. 1-7, majoritatea algoritmilor filtrează un număr de puncte apropiat de cel corect. Excepție fac algoritmii BCAL și Groundfilter, care omit o parte semnificativă din punctele aflate la sol, în toate situațiile. De exemplu, pentru suprafața de probă nr. 4, algoritmul BCAL identifică doar 10 puncte aflate la sol, reprezentând 4 procente din numărul real. Algoritmul Groundfilter are o performanță similară din acest punct de vedere, identificând între 33 și 67 la sută din punctele aflate la sol, în cazul suprafețelor de probă nr. 1-7. Cei doi algoritmi (BCAL și Groundfilter) au totodată cea mai slabă acuratețe de filtrare, din punct de vedere al erorilor medii pătratice. Performanța relativ scăzută a algoritmului Groundfilter în situația terenurilor împădurite confirmă o serie de cercetări anterioare (Kraus și Pfeifer, 1998; Montealegre, 2015).

Tabelul 10.1. Erori ale Modelului Digital Altimetric datorate procesului de filtrare.

Tabelul 10.2. Numărul de puncte filtrate de către algoritmi.

10.1.3. Distribuția spațială a erorilor de filtrare

Calcularea erorilor de filtrare s-a făcut după procedura detaliată în secțiunea 9.4.3, care presupune realizarea diferenței algebrice dintre valorile celulelor a 2 modele raster (cu rezoluția de 0.5 metri):

Modelul Digital Altimetric, generat pe baza rezultatul filtrării și

un Model Digital Altimetric, considerat de referință, obținut pe baza punctelor filtrate corect (pentru suprafața de probă A) sau pe baza măsurătorilor geo-topografice (pentru suprafețele de probă nr. 1-8).

Rezultă prin diferență un model raster al erorilor, suprapus cu celelalte două, unde valoarea fiecărei celule reprezintă eroarea verticală datorată filtrării. Valorile pot fi pozitive sau negative, semnificând supra-estimarea, respectiv sub-estimarea altitudinii de referință. Întrucât interesează mărimea acestor erori, nu neapărat semnul, pentru analiza distribuției erorile au fost transformate în valori absolute.

Erorile absolute au fost apoi stratificate pe patru clase, pentru fiecare clasă calculându-se acoperirea (cu alte cuvinte, ce procent din numărul total de celule al modelului corespunde fiecărei clase). Analiza s-a realizat pentru trei dintre algoritmii de filtrare, câte unul din fiecare categorie: Lasground-new (din categoria adaptare-TIN), MLS (algoritm morfologic) și gLiDAR (algoritm bazat pe interpolare). Analiza distribuției s-a realizat pentru suprafața de probă A, deoarece numărul scăzut de celule ale modelelor pentru suprafețele de probă 1-8 (datorită întinderii mai reduse a acestora) nu garantează semnificația statistică rezultatelor. Limitele stabilite pentru clase, împreună cu rezultatul clasificării, sunt prezentate în tabelul 10.3.

Se observă că suprafața de probă A este acoperită în proporție de 80-85 la sută de erori pe care le putem considera nesemnificative pentru practică (mai mici de 0.20 metri), pentru fiecare dintre algoritmi. Totuși, indiferent de metoda de filtrare, apar și erori majore, de peste 1 metru. Aceste erori acoperă 2-3 procente din suprafață în cazul algoritmilor Lasground-new sau MLS. Performanța algoritmului gLiDAR este relativ slabă, cu aproape 6 procente din suprafață caracterizată prin erori de peste 1 metru.Astfel de erori afectează calitatea modelului suprafeței, fiind necesară corectarea lor anterior folosirii modelului în scopuri practice.

Distribuția spațială a erorilor de filtrare clasificate s-a reprezentat și grafic, sub formă de hărți tematice. În figura 10.2 se prezintă distribuția spațială pentru cazul algoritmului Lasground-new, cu mențiunea că hărțile pentru ceilalți algoritmi sunt similare.

Tabelul 10.3. Clasificarea erorilor verticale de filtrare

Fig. 10.2. Distribuția spațială a erorilor verticale de filtrare pentru algoritmul Lasground-new (de adaugat legenda).

Pentru analiza cantitativă a distribuției spațiale a erorilor s-a apelat la tehnica LISA (Local Indicator of Spatial Association) detaliată în secțiunea 9.4.6. Tehnica urmărește cuantificarea gradului în care valorile tind să se grupeze spațial (prin formarea așa-numitor clustere – grupuri de observații învecinate care au valori similare).

În cazul suprafeței de probă A, analiza LISA scoate în evidență un număr redus de clustere pozitive sau negative (figura 10.3). Majoritatea cluster-elor negative (grupuri de celule învecinate cu valori relativ mici, raportat la media aritmetică a întregii suprafețe) sunt aliniate pe latura nordică a suprafeței de probă A, indicând un posibil efect de margine al tehnicii LISA. Pentru cluster-ele pozitive (grupuri de celule învecinate cu valori relativ mari, raportat la media aritmetică a întregii suprafețe) se observă o anumită tendință de grupare în zona limitrofă albiei râului (care străbate central suprafața de probă A, de la Vest la Est). Această zonă este caracterizată de o versanți abrupți pe ambele părți ale râului și de relieful plat caracteristic albiei (cu alte cuvinte, o discontinuitate topografică – figura 10.4). Este de așteptat ca aceste discontinuități să influențeze procesul de filtrare a norului de puncte, rezultând erori de filtrare relativ ridicate.

Cei trei algoritmi de filtrare analizați au valori similare pentru indicatorul Local Moran’s I (calculat cu formula 9.2), care cuantifică gradul de auto-corelație spațială: 0.28 pentru algoritmul Lasground-new, 0.31 pentru MLS, respectiv 0.32 pentru gLiDAR. Fiind mai apropiate de 0 decât de 1, valorile indică faptul că, la nivelul întregii suprafețe, gradul de organizare pe clustere a erorilor de filtrare este redus.

Fig. 10.3. Distribuția spațială a cluster-elor identificate de tehnica LISA, pentru erorile verticale datorate procesului de filtrare.

Fig. 10.4. Exemplu de discontinuitate a suprafeței terestre, la intersecția dintre zona plată a albiei și versanții abrupți de pe cele două maluri.

10.1.4. Factori de influență a erorilor de filtrare – suprafața de probă A

La analiza procesului de filtrare, s-au avut în vedere o serie de factori cu potențială influență asupra variației erorilor de filtrare (detaliați în secțiunea 9.4.6). Valorile celor cinci factori (înclinarea terenului, Indicele de Rugozitate al Terenului, curbura maximă, gradul de acoperire cu vegetație și densitatea punctelor LiDAR) au fost stratificate pe clase, în funcție de distribuția statistică a valorilor. Semnificația clasificării factorilor și limitele claselor se regăsesc în tabelul 9.4).

Pentru trei dintre algoritmii de filtrare (Lasground-new, MLS și gLiDAR) s-au determinat mediile aritmetice ale erorilor absolute de filtrare, corespunzătoare fiecărei clase a factorilor.

În ceea ce privește înclinarea terenului (figura 10.5) se remarcă o creștere continuă a mediei erorilor absolute de filtrare odată cu înclinarea, pentru fiecare dintre cei trei algoritmi. Algoritmul Lasground-new are o performanță constantă, având cea mai mică valoare a mediei erorilor pentru fiecare dintre categoriile de înclinare. În schimb, algoritmul MLS este aparent cel mai afectat de creșterea înclinării terenului. Performanța acestuia scade relativ la ceilalți doi algoritmi în cazul categoriilor superioare de înclinare (clasele 5 și 6). Valorile indicelui de corelație R2 confirmă această legătură sugerată de reprezentarea grafică (tabelul 10.4).

Situația este asemănătoare pentru Indicele de Rugozitate a Terenului TRI (figura 10.6), precum și pentru curbura maximă a suprafeței (figura 10.7). Performanța fiecărui algoritm de filtrare este afectată negativ de creșterea valorilor celor doi indicatori, care exprimă gradul de fragmentare a suprafeței terenului. Algoritmul Lasground-new are asociate cele mai mici valori pentru media erorilor absolute, la fiecare categorie a celor doi indici. Performanța algoritmului MLS este în scădere relativ la ceilalți doi algoritmi în cazul categoriilor superioare de rugozitate sau curbură maximă a terenului (clasele 4-5). Valorile indicelui de corelație R2 confirmă influența indicatorilor asupra rezultatului filtrării. (tabelul 10.4).

Cu alte cuvinte, se poate spune că acuratețea procesului de filtrare este puternic influențată de condițiile de teren (înclinare, rugozitate, curbură). Din acest punct de vedere, rezultatele confirmă cercetările anterioare (Hyyppa ș.a., 2000; Fisher și Tate, 2006; Meng ș.a., 2010).

Fig. 10.5. Variația mediei erorilor absolute de filtrare pentru cele șase categorii de înclinare a terenului. Liniile verticale delimitează intervalul de încredere al mediei aritmetice.

Fig. 10.6. Variația mediei erorilor absolute de filtrare pentru cele șase categorii ale indicelui TRI. Liniile verticale delimitează intervalul de încredere al mediei aritmetice.

Fig. 10.7. Variația mediei erorilor absolute de filtrare pentru cele șase categorii ale curburii maxime. Liniile verticale delimitează intervalul de încredere al mediei aritmetice.

Creșterea gradului de acoperire cu vegetație al suprafeței (care limitează capacitatea de penetrare la sol a undelor laser) se asociază cu o creștere a mediei erorilor absolute (figura 10.8). În situația terenului descoperit (prima clasă a gradului de acoperire), algoritmii Lasground-new și MLS au o performanță similară. Pentru restul claselor, unde apare vegetația forestieră, se respectă aceeași ierarhie relativă: cele mai mici erori absolute sunt asociate algoritmului Lasground-new, urmat de algoritmii MLS și gLiDAR. Indicele R2, care exprimă legătura dintre erorile de filtrare și gradul de acoperire cu vegetație, este prezentat în tabelul 10.4.

În ceea ce privește densitatea punctelor LiDAR, analiza vizuală a graficului de tip scatterplot nu indică o legătură aparentă între acest factor și valorile erorilor absolute (figura 10.9). Se observă de asemenea un interval de încredere al mediei aritmetice mai semnificativ, în special pentru prima clasă de densitate, datorat numărului redus de observații atribuit acestei clase. Valorile indicelui de corelație R2 indică de asemenea o legătură slabă între densitatea punctelor LiDAR și erorile de filtrare asociate celor trei algoritmi (tabelul 10.4). Sithole și Vosselman (2004) remarcă faptul că influența densității punctelor LiDAR este minoră prin comparație cu efectul complexității suprafeței terestre, iar rezultatele de față confirmă această idee. Și alte cercetări anterioare descoperă o influență relativ scăzută a densității punctelor LiDAR asupra erorilor de filtrare (Zhang și Whitman, 2005; Montealegre ș.a., 2015).

Tabelul 10.4. Valorile coeficientului de corelație R2 pentru legătura între factorii analizați și acuratețea de filtrare a algoritmilor.

Fig. 10.8. Variația mediei erorilor absolute de filtrare pentru cele șase categorii

ale gradului de acoperire cu vegetație.

Liniile verticale delimitează intervalul de încredere al mediei aritmetice.

Fig. 10.9. Variația mediei erorilor absolute de filtrare pentru cele șase categorii

de densitate a punctelor LiDAR.

Liniile verticale delimitează intervalul de încredere al mediei aritmetice.

10.1.5. Factori de influență a erorilor de filtrare – suprafețele de probă nr. 1-7

În ceea ce privește suprafețele de probă nr 1-7, pentru a investiga influența factorilor analizați (înclinarea terenului, Indicele de Rugozitate al Terenului, curbura maximă, gradul de acoperire cu vegetație și densitatea punctelor LiDAR) asupra acurateții de filtrare a algoritmilor s-au calculat:

media aritmetică a fiecărui factor pe suprafață de probă;

eroarea medie pătratică (RMSE) a valorilor de altitudine, pentru fiecare suprafață de probă;

Variația erorii medii pătratice (RMSE) în raport cu fiecare dintre factori s-a reprezentat sub forma graficelor de tip scatterplot. Analiza indică o situație similară cu suprafața de probă A pentru: înclinarea terenului (figura 10.10), Indicele de Rugozitate a Terenului (figura 10.11) și curbura maximă (figura 10.12). În general, suprafețelor care au o medie aritmetică mai mare pentru acești factori li se asociază o scădere a acurateții de filtrare (valori mai mari pentru RMSE).

Referitor la gradul de acoperire cu vegetație, nu este aparentă o legătură între variația acestui factor și erorile medii pătratice asociate suprafețelor (figura 10.13). Această situație este în contrast cu cea întâlnită pentru suprafața de probă A (figura 10.8). O posibilitate este ca variația gradului de acoperire cu vegetație să nu fie suficient de mare pentru a evidenția efectul acestui factor asupra acurateții de filtrare. În fapt, cele șapte suprafețe au un grad de acoperire cu vegetație ridicat (între 78 și 97 la sută).

Efectul modificării densității punctelor LiDAR (figura 10.14) este asemănător cu situația pentru suprafața de probă A. Astfel, între acest factor și eroarea medie pătratică se identifică o corelație slabă.

Valorile indicelui de corelație R2, calculat pentru suprafețele de probă nr. 1-7, sunt prezentate în tabelul 10.4.

Fig. 10.10. Variația erorii medii pătratice a altitudinii pentru suprafețele de probă 1-7,

în raport cu înclinarea medie a terenului.

Fig. 10.11. Variația erorii medii pătratice a altitudinii pentru suprafețele de probă 1-7,

în raport cu media indicelui TRI.

Fig. 10.12. Variația erorii medii pătratice a altitudinii pentru suprafețele de probă 1-7,

în raport cu media curburii maxime a terenului.

Fig. 10.13. Variația erorii medii pătratice a altitudinii pentru suprafețele de probă 1-7,

în raport cu media gradului de acoperire cu vegetație.

Fig. 10.14. Variația erorii medii pătratice a altitudinii pentru suprafețele de probă 1-7, în raport cu media densității punctelor LiDAR.

10.2. Interpolarea suprafeței terenului

10.2.1. Precizări generale

Procesul de interpolare a fost identificat ca una dintre sursele de erori care influențează calitatea Modelului Digital Altimetric. Interpolarea a fost analizată prin testarea a nouă algoritmi de interpolare, pentru suprafața de probă Vidra (figura 9.5). Pentru o parte dintre algoritmi, s-au testat diferite valori ale parametrilor de interpolare (tabelul 4.5). Prin analiza erorilor s-au identificat acele valori care au condus la cea mai precisă reprezentare a suprafeței terestre, adică cea mai mică valoare pentru eroarea medie pătratică.

Acești parametri sunt:

pentru algoritmul Inverse Distance Weighted (IDW): coeficientul u, reprezentând puterea asociată funcției de ponderare (formula 9.8) – cu cât valoarea este mai mare, scăderea influenței punctelor cunoscute odată cu distanța este mai accentuată; în cazul de față valoarea optimă a fost u = 3;

pentru algoritmul Multi-Level B-Spline (BS): parametrul număr maxim de nivele controlează numărul maxim de iterații prin care se optimizează funcțiile B-spline care modelează suprafața; cu cât numărul este mai mare, se generează mai multe suprafețe intermediare, iar teoretic erorile de interpolare scad; în cazul de față, valoarea optimă a fost de 14 nivele, de altfel maximul acceptat de algoritm pentru acest parametru;

pentru algoritmul Cubic Spline (CS): parametrul toleranță controlează gradul de tensiune al suprafeței modelate – cu alte cuvinte, cu cât valoarea scade se reduc schimbările bruște de altitudine ale suprafeței; în acest caz, valoarea finală pentru toleranță a fost de 80;

pentru algoritmul Thin Plate Spline (TPS): nr. maxim de puncte = 20; algoritmul generează o funcție de tip spline pentru fiecare punct nou, folosind punctele cunoscute situate în apropierea acestuia (valoarea standard pentru distanță este de 10 metri). Dacă numărul punctelor apropiate depășește o valoare maximă (în cazul de față 20), se ignoră punctele adiționale, pentru a reduce timpul de procesare;

Valorile de mai sus au condus la interpolarea cu cel mai ridicat grad de acuratețe, pentru fiecare dintre rezoluțiile spațiale analizate (0.5, 1.0, 1.5 și 2.0 metri).

Indicatorii statistici ai erorilor verticale de interpolare sunt prezentați în tabelul 10.5. În toate cazurile, indiferent de algoritm sau rezoluție, erorile medii au valori sub-centimetrice indicând lipsa unei tendințe de supra/sub-estimare a variabilei interpolate. Acest rezultat este în contradicție cu o serie de cercetări anterioare, care descoperă o tendință de supra-estimare a altitudinii în situația terenurilor acoperite cu vegetație forestieră. Astfel, Reutebuch (2003) comunică o eroare medie a Modelului Digital Altimetric de 0.31 metri pentru păduri de conifere, iar Su și Bork (2006) calculează o valoare de 0.23 metri pentru situația pădurilor de plop.

În ceea ce privește eroarea medie pătratică (RMSE), valorile sunt relativ apropiate între algoritmi (de la 0.11 până la 0.40 metri). Ca un alt indicator al acurateții de interpolare, s-a calculat media aritmetică a erorilor de interpolare ca valori absolute (în modul). Pentru o anumită rezoluție, acest indicator are valori apropiate între algoritmi. De exemplu, la rezoluția maximă (0.5 metri) media erorilor absolute are valori între 0.08 și 0.20 metri. În cazul rezoluției minime (2.0 metri), același indicator are valori între 0.26 și 0.30 metri.

Intervalul de valori al erorilor (sau diferența dintre eroarea maximă și eroarea minimă) are valori între 5.88 și 11.29 metri. În general, algoritmilor din categoria Radial Basis Functions (Multilevel B-Spline, Cubic Spline, Thin-Plate Spline și Thin-Plate Spline – TIN) li se asociază un interval al erorilor mai redus. Aceasta deoarece algoritmii de tip RBF urmăresc generarea unei suprafețe netede, caracterizată de o tensiune minimă. Astfel, la nivel local fluctuațiile variabilei interpolate tind să fie reduse.

Tabelul 10.5. Statistici ale erorilor verticale de interpolare

În ansamblu, cel mai puțin precis model al suprafeței rezultă în urma interpolării cu algoritmul Nearest Neighbour (cu valori ale RMSE între 0.28 și 0.41 metri, pentru cele patru rezoluții), urmat de Ordinary Kriging (RMSE între 0.11 și 0.40) și Inverse Distance Weighted (RMSE între 0.20 și 0.36). La polul opus se situează algoritmul Natural Neighbour (NN), suprafețele interpolate de acesta având un grad ridicat de acuratețe (indicatorul RMSE are valori între 0.11 și 0.34 metri). Se face mențiunea că diferențele între algoritmii TPS și TPSTIN sunt foarte reduse, cel puțin din prisma indicatorilor acurateții globale. Având în vedere timpul de procesare mult mai ridicat al TPSTIN (care necesită generarea unui model de tip Rețea Neregulată de Triunghiuri cu milioane de vârfuri), utilizarea acestui algoritm nu se justifică.

Figura 10.15 include reprezentarea grafică a Modelului Digital Altimetric obținut cu algoritmul Nearest Neighbour, la rezoluția de 0.5 metri. Menționăm că, la scara de reprezentare necesară încadrării în pagină a imaginilor, nu se observă diferențe între suprafețele interpolate de diferiți algoritmi. Pentru a remarca diferențele este necesară o scară de reprezentare mai mare, astfel că se prezintă un subset al Modelului Digital Altimetric pentru suprafața de probă Vidra, așa cum a fost interpolat de cei nouă algoritmi (Anexa 5).

Analiza vizuală a rezultatelor evidențiază faptul că fiecare dintre suprafețe include anumite „artefacte” de interpolare (adică modificări sub aspect vizual ale modelului provocate de procedeul de lucru al algoritmilor), în special în zonele unde densitatea punctelor cunoscute este relativ redusă.

Se observă artefacte de interpolare pentru:

algoritmul Nearest Neighbour, unde se distinge forma poligoanelor Thiessen generate anterior interpolării, în zonele unde cresc distanțele dintre puncte cunoscute;

algoritmul Delauney Triangulation, unde se observă triunghiurile rețelei TIN, în zonele cu densitate redusă a punctelor cunoscute;

algoritmul Cubic Spline, caz în care apar artefacte de formă regulată în zonele cu o variație ridicată a altitudinii;

algoritmul Inverse Distance Weighted, care prezintă un efect vizual general de „granularitate”;

algoritmul Ordinary Kriging, unde apare un efect similar cu cel de la IDW, dar mai puțin evident.

În general, se constată ca algoritmul Natural Neighbour, împreună cu algoritmii din categoria RBF (sau spline), generează modele ale suprafeței relativ apropiate de realitate, sub aspect vizual. Acuratețea ridicată de interpolare a algoritmului Nearest Neighbour în situația terenurilor împădurite este constatată și de Sterenczak ș.a. (2016). Performanța relativ scăzută a algoritmului Inverse Distance Weighted este de asemenea în acord cu cercetări anterioare (Bater și Coops, 2009; Ismail, 2016).

Fig. 10.15 Model Digital Altimetric pentru suprafața de probă Vidra, interpolat folosind algoritmul Natural Neighbour (rezoluția modelului: 0.5 metri).

10.2.2. Efectul rezoluției spațiale a modelului asupra acurateții de interpolare

Pentru fiecare dintre algoritmi, erorile verticale datorate interpolării cresc odată cu scăderea rezoluției (adică modificarea în sens pozitiv a dimensiunii celulelor). Erorile medii absolute pentru cele patru rezoluții studiate sunt prezentate în figura 10.16. Se observă că acuratețea de predicție scade la rezoluțiile inferioare. Pentru o anumită rezoluție, ierarhizarea algoritmilor în funcție de acuratețe rămâne în general aceeași. Această influență a rezoluției de modelare este constatată și în cercetările anterioare (Weng, 2006; Bater și Coops, 2009; Ismail, 2016).

Fig. 10.16. Eroarea medie absolută pentru cele patru rezoluții testate.

10.2.3. Distribuția spațială a erorilor de interpolare

Pentru analizarea modului în care erorile de interpolare se distribuie spațial, acestea au fost transformate în valori absolute și grupate pe clase. Limitele claselor de erori și justificarea alegerii acestora este similară cu situația erorilor de filtrare, detaliată în secțiunea 9.4.6.

Rezultatele prezentate în tabelul 10.6 sunt referitoare pentru rezoluția de 0.5 metri, cu mențiunea că situația pentru celelalte trei rezoluții de modelare (1.0, 1.5 și respectiv 2.0 metri) este una foarte asemănătoare. Se remarcă faptul că pentru majoritatea algoritmilor, erorile de interpolare de peste 0.50 metri au o rată redusă de apariție (mai puțin 1 la sută din punctele de validare au asociate astfel de erori). Excepție fac algoritmii Inverse Distance Weighted (erorile mai mari de 0.50 metri sunt asociate la aproximativ 2.5 procente din puncte), Natural Neighbour și Ordinary Kriging (care au erori de interpolare de peste 0.50 metri pentru 6-7 la sută din punctele de validare).

Algoritmii cu acuratețe de interpolare relativ ridicată (Natural Neighbour, Delauney Triangulation, Multi-level B-spline, Thin-Plate Spline, Thin-Plate Spline TIN) au erori de interpolare reduse (de până la 0.20 metri) pentru majoritatea punctelor de validare (93-95 de procente).

Pentru analiza vizuală a distribuției spațiale s-au realizat hărți tematice ale erorilor de interpolare clasificate. Pentru lucrarea de față s-au ales doi algoritmi ca fiind sugestivi: Nearest Neighbour ca algoritm cu acuratețe ridicată de interpolare (figura 10.17), respectiv Ordinary Kriging, care are o performanță scăzută prin raportarea la ceilalți algoritmi de interpolare (figura 10.18).

Tabelul 10.6. Proporția claselor de erori de interpolare pentru suprafața de probă Vidra

Pentru cuantificarea gradului de auto-corelare spațială (măsura în care erori de mărime similară tind să se grupeze) s-a folosit analiza LISA (Local Indicator of Spatial Association). Această tehnică presupune calcularea indicatorului Local Moran’s I, prin care se estimează la nivel global tendința de grupare a valorilor setului de date. Valorile indicatorului Local Moran’s I (formula 9.2) sunt foarte similare între cei nouă algoritmi de interpolare analizați, toate fiind cuprinse în intervalul 0.13-0.17 (Tabelul 10.7). Faptul că valorile sunt apropiate de zero indică un grad redus de organizare spațială pe grupuri (sau clustere) a erorilor de interpolare, indiferent de algoritm. Ca reprezentare grafică a acestei organizări spațiale se prezintă harta cluster-elor pentru algoritmul Inverse Distance Weighted (figura 10.19), cu mențiunea că situația pentru ceilalți algoritmi este foarte similară sub aspect vizual. În reprezentare se observă așezarea aleatoarie a cluster-elor pozitive (grupuri de erori relativ ridicate învecinate între ele), în timp ce majoritatea cluster-elor negative (grupuri de erori relativ scăzute învecinate între ele) se situează:

– în zona sudică a suprafeței de probă Vidra (caracterizată de o înclinare mai redusă și vegetație forestieră mai puțin consistentă sau chiar absentă);

– de-a lungul căii de acces care străbate valea de la Sud la Nord, unde terenul este desigur descoperit;

Tabelul 10.7. Valorile indicatorului Local Moran’s I, pentru erorile de interpolare.

Fig. 10.17. Distribuția spațială a erorilor de interpolare clasificate,

pentru algoritmul Nearest Neighbour.

Fig. 10.18. Distribuția spațială a erorilor de interpolare clasificate,

pentru algoritmul Ordinary Kriging.

Fig.10.19. Distribuția spațială a grupurilor (clustere) de erori identificate prin analiza LISA,

pentru erorile de interpolare.

10.2.4. Factori de influență a erorilor de interpolare

La analiza erorilor de interpolare au fost luați în considerare un număr de cinci factori pentru care s-a analizat o eventuală influență asupra acurateții de interpolare: înclinarea terenului, Indicele de Rugozitate al Suprafeței (TRI), curbura suprafeței, gradul de acoperire cu vegetație și densitatea punctelor cunoscute. Semnificația acestor factori, împreună cu modalitatea de determinare se regăsește în secțiunea 9.5.7. După determinare, factorii au fost clasificați, în funcție de distribuția statistică a valorilor. Limitele claselor și semnificația atribuită acestora sunt prezentate în tabelul 9.6.

Legătura dintre erori și variația factorilor a fost analizată pentru cele patru rezoluții spațiale: 0.5, 1.0, 1.5 și 2.0 metri. Rezultatele fiind foarte asemănătoare, se prezintă doar situația pentru rezoluția de 0.5 metri.

Înclinarea terenului

Analiza erorilor de interpolare, stratificate pe categorii de înclinare a terenului, arată că acuratețea de modelare a suprafeței terenului scade odată cu creșterea înclinării (figura 10.20). De exemplu, pentru algoritmul Natural Neighbour, eroarea medie pătratică (RMSE) crește de la 0.06 metri (pentru punctele de validare aflate pe teren plat) până la 0.51 metri (pentru punctele de validare asociate clasei de înclinare de peste 50 grade). Pentru aceleași două clase, eroarea maximă crește de la 1.38 metri la 4.94 metri. Situația este asemănătoare pentru fiecare dintre cei nouă algoritmi de interpolare analizați. Valorile indicelui de corelație R2 confirmă legătura dintre înclinarea terenului și amplitudinea erorilor de interpolare (tabelul 10.8) pentru fiecare dintre algoritmi. Cea mai strânsă legătură (valori ale indicelui R2 de peste 0.90) se regăsește la algoritmii Nearest Neighbour, Inverse Distance Weighted și Ordinary Kriging. Aceiași algoritmi au cea mai scăzută performanță în ansamblu (tabelul 10.5). Legătura dintre înclinarea terenului și acuratețea procesului de interpolare este în acord cu cercetările anterioare (Hyyppa ș.a., 2005; Aguilar, 2005; Su și Bork, 2006; Ismail, 2016; Sterenczak et al, 2016).

Fig. 10.20. Eroarea medie absolută în funcție de înclinarea terenului.

Pentru claritate, punctele sunt deplasate pe orizontală.

Indicele de Rugozitate a Terenului

Efectul Indicelui de Rugozitate a Terenului (TRI) este asemănător cu cel al înclinării terenului (figura 10.21). De exemplu, pentru indicatorul Ordinary Kriging, valorile erorii medii pătratice cresc de la 0.13 (pentru clasa Fără rugozitate) la 0.80 metri (pentru clasa de Rugozitate ridicată). Valorile indicelui de corelație R2 (peste 0.85 pentru fiecare algoritm) confirmă efectul indicelui TRI asupra variației erorilor de interpolare (tabelul 10.8). La fel ca în cazul înclinării terenului, cea mai strânsă legătură apare pentru cei trei algoritmi cu cea mai slabă acuratețe de interpolare (NeN, IDW și OK), cu valori ale indicelui R2 de 0.97, 0.95 și respectiv 0.94.

Curbura maximă a suprafeței

La fel ca înclinarea terenului sau Indicele de Rugozitate, curbura maximă exprimă într-un anumit sens gradul de complexitate al suprafeței terestre. Astfel, la fel ca în cazul primilor doi indicatori, erorile de interpolare sunt în general mai mari pentru clasele de curbură ridicată (figura 10.22). Modificarea cea mai accentuată apare între ultimele două clase ale curburii. Ca exemplu, eroarea medie pătratică (RMSE) pentru algoritmul Thin-Plate Spline are valori de 0.12-0.19 metri pentru primele patru clase de curbură, în timp ce valoarea RMSE pentru cea de-a cincea clasă (curbură foarte ridicată) este de 0.47 metri. Indicele de corelație R2 are valori de peste 0.85 pentru fiecare algoritm, indicând o legătură puternică între curbura suprafeței și erorile de interpolare.

Fig. 10.21. Eroarea medie absolută în funcție de Indicele de Rugozitate a Terenului (TRI).

Pentru claritate, punctele sunt deplasate pe orizontală.

Fig. 10.22. Eroarea medie absolută în funcție de curbura maximă a terenului.

Pentru claritate, punctele sunt deplasate pe orizontală.

Gradul de acoperire cu vegetație

Vegetația limitează capacitatea de penetrare la nivelul terenului a undelor laser, afectând nu doar densitatea punctelor aflate la sol ci și modul de distribuție al acestora (densitatea va fi ridicată în cazul terenului descoperit și mult mai redusă în zonele acoperite de pădure). Din acest motiv este de așteptat ca acest factor să influențeze într-o anumită măsură variația erorilor de interpolare. Rezultatele confirmă în parte această ipoteză (figura 10.23). Pentru o parte dintre algoritmi (Nearest Neighbour, Ordinary Kriging, Inverse Distance Weighted) erorile cresc pentru clasele superioare ale gradului de acoperire cu vegetație. Însă ceilalți algoritmi au o performanță relativ stabilă din acest punct de vedere. De exemplu, valorile erorii medii pătratice (RMSE) pentru algoritmul Delauney Triangulation variază între 0.10 și 0.16 metri. Pentru algoritmul Cubic Spline se remarcă o oarecare scădere a RMSE (de la 0.22 metri pentru clasa Teren descoperit la 0.15 metri pentru clasa de Acoperire foarte ridicată). Acest comportament este însă cel mai probabil aleatoriu, fără semnificație statistică. Valorile indicelui de corelație R2 sunt prezentate în tabelul 10.8.

Fig. 10.23. Eroarea medie absolută în funcție de gradul de acoperire cu vegetație.

Pentru claritate, punctele sunt deplasate pe orizontală.

Densitatea punctelor de predicție

În ceea ce privește densitatea punctelor din setul de date pentru predicție (folosit la interpolarea suprafețelor), se remarcă o legătură puternică între acest factor și erorile de interpolare (figura 10.24). Influența este mai puternică pentru algoritmii care au în general o acuratețe de interpolare scăzută (Inverse Distance Weighted, Nearest Neighbour, Ordinary Kriging). Ceilalți algoritmi au valori relativ apropiate ale erorii medii pătratice (RMSE) indiferent de clasa de densitate a punctelor. De exemplu, indicatorul RMSE pentru algoritmul Natural Neighbour are o variație redusă, de la 0.14 la 0.11 metri. Valorile indicelui de corelație R2 (tabelul 10.8), indică o corelație negativă pentru fiecare algoritm: creșterea densității punctelor folosite la predicție se asociază unor valori reduse pentru eroarea medie pătratică. Efectul densității punctelor cunoscute asupra preciziei de interpolare confirmă cercetările anterioare (Aguilar, 2005; Guo, 2010). De asemenea se observă o mai mare variabilitate între acuratețea algoritmilor, în situația claselor de densitate scăzută. Acest aspect este remarcat și de Chaplot (2006).

Tabelul 10.8. Valoarea coeficienților de corelație R2 pentru legătura între factorii analizați și acuratețea de filtrare a algoritmilor.

Fig. 10.24. Eroarea medie absolută în funcție de densitatea punctelor cunoscute.

Pentru claritate, punctele sunt deplasate pe orizontală.

11. CONCLUZII, CONTRIBUȚII ȘI RECOMANDĂRI

11.1. Concluzii

Teza de doctorat, intitulată „…” a avut ca obiect general stabilirea potențialului tehnologia de scanare laser (LiDAR) ca sursă de date despre suprafața terestră, în condițiile terenurilor acoperite cu vegetație forestieră.

Prin parcurgerea etapelor de cercetare necesare atingerii obiectivului propus, s-au desprins o serie de remarci generale, organizate pe secțiuni:

Stadiul actual al cunoștințelor

Metodele moderne ale geomaticii forestiere permit o cunoaștere la nivel de detaliu a pădurii, ca poziție, structură, caracteristici și evoluție în timp.

Tehnologia LiDAR, bazată pe emisia undelor electromagnetice de tip laser, se înscrie în metodele moderne de cartare ale geomaticii, alături de fotogrametrie sau teledetecție.

Pe scurt, tehnologia folosește fascicule coerente de lumină pentru a determina indirect distanța între emițător și suprafața care produce reflexia fasciculelor, cu scopul de a poziționa în spațiu suprafața scanată.

Avantajul tehnologiei LiDAR pentru silvicultură este dat de capacitatea de penetrare a frunzișului, care permite pulsurilor laser să atingă suprafața terestră. Se pot astfel obține simultan date despre caracteristicile vegetației și ale reliefului.

Scanarea laser aeriană (ALS), care presupune montarea senzorului LiDAR pe o platformă aeropurtată, aduce tehnologia pe același nivel cu fotogrametria, din punct de vedere al randamentului de acoperire a suprafețelor.

Tehnologia LiDAR s-a răspândit pe plan internațional în ultimele decenii, ca o metodă precisă de cartare a suprafețelor la nivel de detaliu. Această răspândire s-a realizat și în domeniul geomaticii forestiere, unde scanarea laser permite automatizarea procesului de obținere a datelor geografice și structurale ale pădurii și solului.

În România, tehnologia este încă în faza de pionerat, studiile realizate până în prezent fiind axate în general asupra estimării anumitor caracteristici dendometrice.

Datele obținute prin scanare sunt reprezentate sub forma unui nor de puncte pentru care se cunoaște poziția planimetrică și altimetrică. Trecerea la o structură comună de reprezentare a suprafeței terestre, sub forma unui Model Digital Altimetric, este un proces complex, afectat de multiple surse de eroare.

Dintre sursele de erori care intervin pe parcursul procesului de trecere de la norul de puncte inițial la un model al suprafeței terestre, s-au ales două pentru analiză: filtrarea datelor LiDAR (prin care se separă din norul de puncte acele observații aflate la nivelul solului) și interpolarea (prin care se generează o suprafață continuă pe baza observațiilor punctuale).

Caracteristici ale condițiilor de cercetare.

Scopul stabilit pentru teza de doctorat a fost evaluarea potențialului tehnologiei LiDAR ca sursă de date altimetrice pentru reprezentarea cu precizie ridicată a suprafeței terestre, în situația terenurilor împădurite

Pentru atingerea scopului propus, s-au stabilit două areale de studiu: Mălaia și Vidra caracterizate de un relief montan și acoperire extinsă cu vegetație forestieră.

Pentru cercetare s-au folosit două tipuri de date: date LiDAR colectate prin scanare laser aeriană și date geo-topografice, colectate folosind stație totală și echipament GPS.

Metodologia de cercetare.

Pentru estimarea acurateții de modelare a suprafeței terestre, s-au identificat două tipuri de erori (de filtrare și respectiv interpolare) care au fost definite, calculate și analizate stastistic.

Pentru erorile de filtrare, s-a testat un număr de nouă algoritmi de filtrare. S-au analizat condițiile de filtrare pentru două cazuri: unul pentru suprafața de probă Mălaia (unde referința este setul de date LiDAR corect clasificat), cel de-al doilea pentru șapte suprafețe de probă delimitate în teren (unde referința este setul de date independent rezultat prin prelucrarea măsurătorilor topografice).

Factorii identificați ca având o posibilă influență asupra erorilor de filtrare sunt: înclinarea terenului, Indicele de Rugozitate al Suprafeței (TRI), curbura în profil a suprafeței, gradul de acoperire cu vegetație și densitatea punctelor LiDAR.

Pentru erorile de interpolare, s-a testat un număr de nouă algoritmi de interpolare a suprafeței. Procesul de interpolare s-a analizat pentru cazul suprafeței de probă Vidra, iar ca factori cu influență asupra erorilor s-au identificat: înclinarea terenului, gradul de acoperire cu vegetație și densitatea punctelor LiDAR. Analiza interpolării s-a efectuat pentru trei rezoluții spațiale de modelare: 0.5, 1.0 și 1.5 metri.

Rezultate și discuții

Cu privire la filtrarea datelor LiDAR

Dintre algoritmii de filtrare analizați, acuratețea de filtrare cea mai bună este oferită de algoritmul Lasground-new – eroare medie pătratică de 0.34 metri pentru suprafața de probă A și un procent de 84 la sută din suprafață caracterizat de erori de determinare a altitudinii între 0.0-0.20 metri.

Principiul de filtrare al algoritmilor afectează acuratețea filtrării: cei trei algoritmi bazați pe interpolare (GroundFilter, BCAL și MCC) se asociază cu erori medii pătratice mai mari, relativ la algoritmii morfologici sau cei din categoria densificare-TIN.

Șapte dintre cei nouă algoritmi analizați filtrează un număr de puncte mai mare decât cel real – erorile de comitere (atribuirea în clasa Teren a punctelor aflate deasupra suprafeței terestre) sunt deci mai probabile decât cele de omitere.

Distribuția spațială a erorilor a fost analizată pentru trei algoritmi de filtrare, câte unul din fiecare categorie (Lasground-new, MLS, gLiDAR) – rezultatele indică o acuratețe ridicată pentru cea mai mare parte a suprafeței de probă A (între 74 și 84 la sută din suprafață este caracterizată de erori absolute de până la 0.20 metri). Totuși, apar și erori semnificative, de peste 1.00 metri, chiar dacă acestea sunt rare (între 2 și 6 la sută din suprafața de probă, în funcție de algoritm). Se impune astfel necesitatea corectării manuale a rezultatului filtrării pentru a putea obține o reprezentare adecvată a suprafeței terestre.

Situația erorilor pentru suprafețele de probă 1-7, cartate folosind metode geo-topografice, este în linii mari aceeași cu situația descrisă pentru suprafața de probă.

Condițiile de teren (înclinarea terenului, Indicele de Rugozitate a Suprafeței și curbura maximă a suprafeței) au o influență semnificativă asupra erorilor de filtrare (indicele de corelație R2 are valori între 0.88 și 1.00, în funcție de factor și algoritm). Astfel, rezultatele confirmă cercetările anterioare.

Pentru analiza corelației spațiale a erorilor s-a calculat indicatorul Local Moran’s I; valorile acestuia indică un grad redus de grupare pe cluster-e a erorilor de filtrare.

În final, referitor la procesul de filtrare automată a datelor LiDAR, putem afirma că rezultatele sunt mulțumitoare sub aspectul preciziei. Trebuie ținut seama de faptul că scanarea laser s-a realizat în perioada de vegetație, iar ecosistemele forestiere reprezintă o provocare pentru orice algoritm de filtrare (Guan ș.a., 2014; Maguya ș.a., 2014; Montealegre ș.a., 2015).

Cu privire la interpolarea suprafeței terestre

Cei nouă algoritmi de interpolare testați au o precizie similară, estimată prin eroarea medie pătratică. Media aritmetică a erorilor este apropiată de zero în toate cazurile, indicând lipsa unei tendințe de supra-estimare sau sub-estimare a predicției valorilor de altitudine.

Întrucât nivelul de acuratețe globală este relativ similar între algoritmi, se impune analiza calitativă (sub aspect vizual) a suprafețelor interpolate (Wood și Fischer, 1993). Se remarcă prezența artefactelor vizuale de interpolate în toate cazurile, dar acestea sunt mai puțin pronunțate în cazul algoritmului Nearest Neighbour și al algoritmilor din categoria RBF (Radial Basis Functions).

Se constată o modificare semnificativă a erorilor de interpolare în funcția de rezoluția modelului. Trebuie ținut seama însă de faptul că la estimarea erorilor s-a folosit metoda cros-validării – în acest caz, valoarea unui punct de validare se compară cu cel mai apropiat punct pentru care s-a realizat predicția (reprezentat de centrul celulei modelului); atunci când rezoluția scade (celulele au o dimensiune mai mare), distanțele între punctele de validare și cele de predicție cresc, deci și posibilitatea unei diferențe reale (nu cauzate de o eroare) între valorile pentru cele două puncte.

Modificarea înclinației terenului influențează negativ acuratețea de predicție indiferent de algoritmul de interpolare folosit, rezultat aflat în acord cu cercetările anterioare (Aguilar, 2005; Su și Bork, 2006; Sterenczak et al, 2016). Astfel, dacă pentru terenurile plate (cu înclinare de până la 5 grade) eroarea medie absolută este de sub 0.10 metri, aceasta crește la peste 0.50 metri în zonele cu înclinație de peste 30-40 grade.

Densitatea punctelor cunoscute influențează de asemenea acuratețea de interpolare, rezultat care confirmă cercetări anterioare (Aguilar, 2005; Guo, 2010). Se remarcă faptul că, atunci când densitatea punctelor scade, se accentuează diferențele de acuratețe între algoritmi (vezi și Chaplot, 2006). Astfel, atunci când datele folosite la interpolare nu o densitate corespunzătoare, capătă importanță alegerea algoritmului de interpolare.

În ceea ce privește gradul de acoperire cu vegetație, acest factor influențează negativ magnitudinea erorilor. Astfel, la acoperire de până la 20 la sută erorile medii absolute au valori în jur de 0.20 metri, la o acoperire cu vegetație de peste 90 la sută erorile medii absolute au valori de 0.35-0.45 metri.

În ansamblu, luând în considerare valorile erorilor dar și analiza vizuală a rezultatelor, cel mai adecvat algoritm de interpolare pentru situația de față este Natural Neighbour. Dincolo de erori, care sunt similare cu cele ale celorlalți algoritmi, avantajul Natural Neighbour este timpul redus de procesare și ușurința utilizării (nu necesită setarea de parametri).

În final, putem afirma că acuratețea de interpolare este relativ bună, ținând seama de topografia complexă a suprafeței de studiu. Totuși, apar izolat erori de ordinul metrilor, care necesită corecții manuale sau automate, de exemplu aplicarea filtrelor de trecere sau “umplerea golurilor” (sink-fill).

În ansamblu, putem concluziona că procesarea automată a datelor LiDAR cu scopul de a obține un model detaliat al suprafeței terestre oferă un grad satisfăcător de acuratețe. Trebuie avut în vedere faptul că datele LiDAR folosite la cercetare nu au fost înregistrate în scopuri forestiere. Astfel, prin planificarea adecvată scopurilor a campaniei de zbor (de exemplu, realizarea înregistărilor în perioada de toamnă/primăvară, când frunzișul este absent dar terenul nu este acoperit de zăpadă) ne așteptăm ca volumul de muncă implicat la procesare să se reducă, în timp ce acuratețea de reprezentare a reliefului să crească.

Stadiul actual al dezvoltării algoritmilor (atât în cazul procesului de filtrare, cât și al interpolării) impune realizarea unor corecții manuale a rezultatelor. Alegerea algoritmului, stabilirea parametrilor de operare și condițiile de teren influențează volumul de lucru necesar la corecția manuală. Dintre cele două operațiuni, filtrarea este cea care are potențialul de a produce erori verticale mai mari. Astfel, în urma filtrării folosind algoritmul Lasground-new (cel mai precis rezultat) se obține o reprezentare a suprafeței terestre caracterizată în proporție de 84 la sută de erori verticale mai mici de 0.20 metri și mai de 2 procente cu erori majore, de peste 1 metru. În situația interpolării, aplicarea algoritmului Nearest Neighbour (cel mai precis rezultat) conduce la un model al terenului cu erori mai mici de 0.20 metri pentru aproape 94 la sută din suprafață, respectiv erori majore (de peste 1 metru) cu o acoperire insignifiantă (0.02 procente din suprafața analizată).

11.2. Contribuții personale

Rezultatele obținute prin desfășurarea cercetărilor și concluziile formulate pe baza acestora scop în evidență o serie de contribuții personale, cu aspect de noutate în cercetarea realizată în țara noastră:

Realizarea unei cercetări bibliografice comprehensive asupra tehnologiei LiDAR, ca metodă generală a geomaticii dar în special din perspectiva aplicabilității în domeniul forestier.

Axarea cercetărilor, în premieră la noi în țară, pe potențialul scanării laser de a oferi date privind relieful suprafeței terestre, în situația terenurilor acoperite cu vegetație forestieră.

Stabilirea unei proceduri ample de evaluare a erorilor, ca distribuție statistică și spațială, care permite compararea obiectivă a algoritmilor sau procedeelor de lucru. Această metodă a inclus de asemenea realizarea măsurătorilor geo-topografice, ca sursă de date independente.

Testarea detaliată a celor două etape principale de prelucrare a datelor LiDAR (filtrarea, respectiv interpolarea suprafeței): s-au analizat un număr ridicat de algoritmi, cu diferite procedee de lucru.

Analizarea erorilor nu doar sub aspectul distribuției statistice, dar și ca distribuție spațială. S-au stabilit de asemenea proceduri de estimare a influenței condițiilor de teren sau vegetație asupra procesului de filtrare, respectiv interpolare.

Formularea avantajelor tehnologiei LiDAR, dar și identificarea limitărilor, din perspectiva geomaticii forestiere.

Definirea, pe baza experienței obținute, a unor direcții promițătoare pentru cercetării viitoare, care să extindă aplicabilitatea în țara noastră a tehnologiei LiDAR, pentru domeniul forestier.