Geostatistica

Ricostruzione tramite GIS della superficie piezometrica della falda alluvionale nell’Alta valle del Tevere

6. Analisi dei livelli di falda con approccio stocastico

6.2 Geostatistica

La geostatistica è quella branca della statistica che si occupa dell’analisi di dati spaziali [Raspa, 1995].

Il suo campo classico di applicazione sono le Scienze della Terra, in particolar modo nella Geologia, Geologia Ambientale, Ecologia, Meteorologia, Agronomia.

La geostatistica si occupa di valutare l’autocorrelazione spaziale dei dati, cercando di verificare se osservazioni effettuate su punti vicini presentino effettivamente una minore variabilità rispetto ad osservazioni effettuate tra punti distanti. L’obiettivo è quindi valutare l’effetto della posizione del punto di misura sulla variabilità del dato osservato. Tale variabilità viene di solito elaborata con lo strumento del semivariogramma, che valuta il grado di variabilità di punti a distanze crescenti.

Lo studio della variabilità spaziale è necessario per la successiva fase di predizione spaziale con cui si possono fornire delle stime sul valore assunto e sugli errori commessi da una variabile in una posizione in cui la misurazione non è stata effettuata. Il metodo di interpolazione statistico più noto è il kriging.

Esistono vari tipi di kriging:

• Kriging ordinario (per processi stazionari)
• Kriging universale (per processi non stazionari)

Il dottor Herbert Sichel ed il professor Danie Krige sono stati i pionieri della geostatistica, mentre Georges Matheron è unanimamente riconosciuto come il padre della Geostatistica.

6.2.1 Autocorrelazione

L’autocorrelazione definisce il grado di dipendenza spaziale tra i valori assunti da una variabile campionata.

Una caratteristica intuitiva dell’ambiente è che le sue proprietà sono in relazione fra di loro in una qualche scala, grande o piccola che sia. Questa situazione è definita autocorrelazione (Figura 6.3).

Questo significa che valori campionati in luoghi vicini tra di loro, tendono ad avere comportamenti simili, mentre valori di una stessa variabile misurati in campioni raccolti in luoghi lontani tra di loro tendono ad avere comportamenti differenti, o almeno tendono a differire dai valori medi che si riscontrano nei due luoghi stessi. In tal senso, la correlazione fra i valori della variabile tende a diminuire con l’aumentare della distanza.

Adottando il punto di vista stocastico, ogni punto dello spazio non ha un solo valore per una proprietà, ma un intero insieme di valori. In questo modo il valore osservato diventa un valore estratto a caso da un infinito numero di valori possibili, assunti da qualche distribuzione di probabilità per una qualche legge. Questo significa che in ogni punto dello spazio c’è una variazione. Così, ad ogni punto x0, una proprietà, Z(x0), è trattata come una variabile casuale, generalmente continua e non discreta, con una

media , una varianza (), momenti di ordine più alto ed una funzione di probabilità di densità cumulativa. Questa variabile ha quindi una distribuzione di probabilità da cui viene estratto il valore reale. Il set di valori reali (misurati) di Z può considerarsi come una variabile regionalizzata in quanto, è una variabile il cui valore è fortemente condizionato dalla posizione spaziale. Detta Z(x) la variabile regionalizzata nella posizione x, potremmo scrivere:

Z (x) = a+R(x)

in cui a è la componente aleatoria e R(x) la componente regionalizzata. Quando a è dominante rispetto a R(x) si studia la variabile con i metodi classici della statistica. In caso contrario si ricade nel campo della geostatistica.

I valori delle variabili regionalizzate tendono ad essere in relazione tra loro e, come detto, due valori vicini ad un terzo tendono ad essere simili, mentre quelli più lontani lo sono meno.

6.2.2 Ipotesi di stazionarietà

Un processo stazionario è un processo stocastico in cui la funzione di densità di probabilità di una qualche variabile casuale Z non cambia né nel tempo né nello spazio. Quindi anche i parametri media e varianza non cambiano nel tempo e nello spazio, di conseguenza la distribuzione del processo casuale mantiene gli stessi attributi ovunque.

Supponiamo di misurare una variabile regionalizzata Z(x) in una particolare direzione dello spazio ottenendo alcuni valori per ogni punto in questa direzione. Si potrebbe ottenere un andamento di valori simile a quello mostrato.

Come si vede, è presente un cambiamento sistematico del valore medio della variabile, rilevante alla scala dell’osservazione, e allora non si può più assumere una stazionarietà della media, o stazionarietà del modello. Questa condizione di non stazionarietà dei valori della variabile e data dalla somma di due componenti:

• drift (aumento sistematico del valore medio)
• residuo (variabilità attorno alla componente sistematica)

Il drift è essenzialmente il valore medio della variabile in funzione della localizzazione in cui la variabile è stata misurata e può essere scritto come m(x). Il valore della variabile Z per ogni punto x può essere espresso allora come la somma della media del punto x e la deviazione dalla media del valore misurato, Y(x):

Z(x) = Y(x) + m(x)

Le assunzioni fatte sono le seguenti:

• il residuo Y(x) è una variabile casuale stazionaria, in cui la media è costante
• il drift sia una funzione deterministica della localizzazione
• il drift e il residuo sono non correlati

6.2.3 Semivariogramma

Il semivariogramma è un algoritmo geostatistico che viene impiegato per valutare l’autocorrelazione spaziale di dati osservati in punti georeferiti.

La funzione semivariogramma interpola la semivarianza dei valori osservati in gruppi di coppie di punti a determinate distanze (lag) secondo una certa direzione.

La semivarianza è pari a:

in cui Z è il valore di una misura in un particolare punto, h è un intervallo di distanza tra punti di misurazione (lag) e m(h) rappresenta il numero di coppie di osservazioni effettuate alla distanza h.

Gli assi del semivariogramma sperimentale sono distanze tra coppie di dati (asse x) e semivarinanza (asse y) (Figura 6.5). Il grafico ottenuto dall’applicazione della formula 6.1 è formato da una serie di punti distinti, distanti tra loro un certo lag definito. Il semivariogramma sperimentale (punti discreti) deve essere interpolato con diverse funzioni matematiche, in modo da determinare il tipo di autocorrelazione spaziale della variabile misurata, che saranno successivamente usate nell’interpolazione con il kriging. Nel processo di kriging è infatti richiesta una funzione continua per l’assegnamento dei pesi in tutti i punti, che è legata al valore della semivarianza. Tale modello di semivariogramma si desume a partire dall’osservazione del semivariogramma sperimentale.

Si può immaginare una funzione che meglio approssimi i punti del variogramma sperimentale; il problema è però capire quanto questa curva seguirà i dati sperimentali. Per modellare bene il variogramma bisogna quindi porre attenzione sull’andamento generale dei punti e non sulle singole fluttuazioni, quindi il tipo di funzione da scegliere dovrà essere la più semplice possibile, in relazione sempre alla complessità dell’andamento del variogramma.

I parametri stimati sono:

• Nugget: descrive il livello di variabilità casuale;
• Sill: valore massimo della semivarianza quando si ha stazionarietà (esso approssima per eccesso la varianza campionaria);
• Range: rappresenta la distanza massima entro la quale si manifesta correlazione tra semivarianza e lag;

La funzione semivariogramma è una funzione di correlazione spaziale, con determinate caratteristiche grafiche e quantitative, dette proprietà strutturali, quali:

• simmetria
• continuità
• zona d’influenza
• comportamento vicino all’origine
• anisotropia
• drift

La modellazione del variogramma è un passo molto delicato in quanto consiste nello scegliere la curva che meglio approssima la serie discreta di punti ottenuti dall’analisi del variogramma sperimentale. La procedura di adattamento dei modelli (fitting model) al variogramma è difficile per varie ragioni:

1. l’accuratezza delle semivarianze osservate non è costante
2. il variogramma può contenere molte fluttuazioni punto a punto

Perciò è raccomandabile una procedura di fitting del modello che inglobi sia l’aspetto visuale nell’andamento della funzione che quello statistico, come segue:

1. va tracciato il variogramma sperimentale
2. si sceglie uno o più modelli tra quelli disponibili, i quali abbiano approssimativamente una forma corretta e con sufficiente dettaglio, in modo da rispettare le tendenze principali dei valori sperimentali
3. si adatta il modello attraverso procedure statistiche di minimizzazione degli errori
4. si ispeziona il risultato grafico, per valutare qualitativamente la bontà del risultato

Esistono vari tipi di modelli per l’approssimazione dei variogrammi sperimentali, tra i più usati troviamo:

• il modello esponenziale
• il modello sferico
• il modello gaussiano
• il modello lineare

6.2.4 Calcolo del variogramma

La stima della funzione variogramma viene effettuata sulla base dei dati provenienti dal campionamento del fenomeno oggetto di studio. Se si hanno dati campionati secondo una maglia regolare il calcolo è molto semplice poiché, data la stazionarietà dell’incremento Z(x+h) – Z(x), risulta immediato calcolare la funzione variogramma per una certa direzione e per un determinato lag h. Il calcolo del variogramma si basa sul considerare le differenze dei valori della variabile regionalizzata in due localizzazioni differenti, separate da una distanza h. La procedura da seguire è questa:

1. Si parte facendo la differenza tra i valori z(x1) e z(x2), poi tra z(x2) e z(x3), fino alla coppia z(xi-1) e z(xi). Le differenze saranno uguali a m(x), ovvero al numero di coppie di campioni per questo lag
2. Il risultato di ogni differenza si eleva al quadrato
3. Si sommano tutti i quadrati
4. Si divide questa somma per 2m(h), come da formula 6.1
5. Si ripete la stessa procedura da 1 a 4 per il secondo lag (il doppio del primo)
6. Si ripete la stessa procedura da 1 a 4 per il terzo lag (tre volte il primo)
7. Si ripete la procedura fino all’ultima distanza di lag che si vuole

Generalmente si sceglie un numero di lag che varia da 10 a 20. Infatti, un basso numero di lag (minore di 10) potrebbe far perdere di significato al variogramma, abbassando il potere risolutivo del variogramma. Al contrario, un alto numero di lag porterebbe a considerare troppi lag oltre il range massimo di correlazione spaziale tra i dati.

Purtroppo però non sempre si ha una maglia regolare di campionature, bisogna allora aggiustare la tecnica per il calcolo del variogramma. In questo caso infatti potrebbe accadere che in una certa direzione r non cada nessun punto campionato, perché disposti in modo irregolare. Per ovviare a ciò allora si considera una direzione r, individuata tramite un angolo a partire da una direzione di riferimento, con una tolleranza angolare Δ. Una tolleranza Δh deve anche essere data sulla distanza. Così facendo tutte le coppie aventi distanza compresa tra h-Δh e h+Δh e allineate secondo la direzione compresa tra -Δ e +Δ contribuiscono al calcolo del variogramma.

I valori delle tolleranze da adottare dipendono ovviamente dalla quantità di campioni di cui si dispone: più sono numerosi i campioni e più piccole possono essere le tolleranze, consentendo una maggior precisione nel calcolo dei variogrammi.

6.2.5 Mappa del variogramma

La mappa del variogramma è un grafico che rappresenta l’andamento della varianza nello spazio. Viene costruito, utilizzando una griglia con celle quadrate (o con settori circolari), sulla base dei valori delle ordinate del variogramma sperimentale. L’ampiezza delle celle è determinata dall’ampiezza dei lag, più è alto il numero di lag tanto più le celle sono piccole.

La mappa del variogramma fa sì che sia subito comprensibile come vari nello spazio l’anisotropia per poter decidere su quali direzioni principali concentrare l’attenzione con il variogramma. Ogni cella della mappa è rappresentata infatti con una simbologia di colori variabile a seconda del valore del variogramma secondo una scala cromatica propria (Figura 6.8).

Questo grafico, come anticipato, permette di individuare la direzione di maggiore continuità del dato semplicemente osservandolo dal centro ed individuando la direzione per cui si hanno valori più bassi di varianza (quadratini rosa scuro). Questo differente andamento della variabilità spaziale è noto come anisotropia. L’anisotropia è la proprietà per la quale un determinato oggetto ha caratteristiche che dipendono dalla direzione lungo la quale esse sono considerate.

6.2.6 Cross-Validation

La cross-validation è una procedura per cui iterativamente ogni campione viene escluso dal dataset ed interpolato attraverso il valore degli altri, utilizzando il modello di variogramma che si vuole testare. Il confronto tra il valore stimato ed il valore reale è detto residuo della cross-validation.

Lo studio del dataset dei residui ci indica il comportamento del modello. I principali parametri da studiare sono:

1. La media dei residui, che indica l’accuratezza della stima, il quale deve essere prossimo a zero
2. L’errore minimo quadrato (RMSE), che indica la precisione , il quale deve essere il più piccolo possibile
3. La deviazione standard del kriging (MSDR)  , la quale permette di confrontare la grandezza degli errori predetti con gli errori realmente commessi, si confrontano cioè la varianza vera e quella calcolata nella cross-validation.

dove Zi è il valore stimato i-esimo, Z’ il valore reale,  la varianza dell’errore ed n il numero di campioni.

6.2.7 Kriging e Cokriging

Per poter ricostruire la superficie cercata è necessario interpolare i dati disponibili per stimare i valori dove non si hanno campioni. In questo senso il kriging ci viene incontro, in quanto dà una soluzione al problema della stima basato su un modello continuo di variazione spaziale stocastica. Esso fa il miglior uso della conoscenza della variabile, prendendo in considerazione il modo in cui una proprietà varia nello spazio attraverso il modello del variogramma che è stato scelto e validato.

Esistono diversi tipi di kriging, tra cui kriging ordinario (ordinary kriging) e kriging universale (universal kriging), per diverse tipologie di variabili. Ciò che li differenzia è il tipo di variabile usata: il kriging ordinario può lavorare solo con variabili stazionarie del secondo ordine (presentano cioè media costante e varianza dipendente solo dal lag muovendosi da punto a punto), il kriging universale può invece lavorare anche con variabili non stazionarie (che presentano cioè un drift). Come anticipato una delle assunzioni fatte nel kriging ordinario è la stazionarietà del dato da stimare. Questo significa che muovendosi da una zona ad un’altra del campo la media dei valori è pressoché costante. Quando invece esiste un significativo trend spaziale del dato (caratteristica intrinseca del dato stesso che fa si che la media dei valori non sia costante ma vari da punto a punto) questa assunzione viene meno. La condizione di stazionarietà del dato può essere comunque ristabilita attraverso l’introduzione di una funzione deterministica che descriva il drift, cioè l’andamento della media, in modo da poter isolare il residuo, cioè la parte aleatoria del dato. Il kriging universale quindi modella e sottrae il drift presente nel dato tramite una funzione deterministica, ed analizza la sola componente aleatoria (residuo).

Dato che il drift ed il residuo sono non correlati è possibile modellare il drift tramite una funzione deterministica della localizzazione, si può scrivere quindi il valore della variabile Z (x) come:

dove  è la somma di un set di funzioni f k (x) polinomiali di ordine 1 o 2 ed f 0 (x)=1 e il termine m(x) è il componente stocastico (variogramma).

L’universal kriging, utilizzato in condizioni di non stazionarietà, fornisce la stima puntuale (nel punto x 0) attraverso la formula:

la quale fornirà una stima corretta solo se si rispetta la condizione:

La soluzione del sistema qui di seguito permette di quantificare i pesi  per la stima del generico punto:

dove ( xi , x j) è la semivarianza dei residui tra i punti xi e xj mentre (xi , x0)è la semivarianza tra l’i-esimo punto campionato ed il punto bersaglio della stima x0.

Nel kriging, sia ordinario che universale, come appena visto, tramite un sistema di equazioni lineari, si assegnano i pesi per i campioni circostanti al punto analizzato che minimizzano la varianza totale degli errori. Solitamente i quattro cinque campioni più vicini contribuiscono per l’80% del peso totale, il restante 20% viene assegnato all’incirca ad altri dieci punti vicini. Si hanno vari fattori che influenzano i pesi tra cui:

1. i campioni vicini portano più peso di quelli lontani e l’entità del peso dipende dalla loro posizione e dal variogramma;
2. i campioni raggruppati in un certo intorno portano meno peso di quelli isolati.

Nell’ambito della statistica multivariata il metodo di interpolazione usato è il cokriging. Questo metodo è molto simile al ordinary kriging, l’unica differenza risiede nel fatto che, nello stimare la variabile principale (target) non si basa solo sui valori della variabile esaminata ma prende in considerazione anche altre variabili ausiliarie.

Si consideri il caso in cui si abbiano due variabili spazialmente correlate Z1(x) e Z2(x), campionate rispettivamente in un set di n1 e n2 locazioni, volendo effettuare la stima del valore della variabile Z1(x0) tramite i valori circostanti di Z1(x) e Z2(x) utilizzeremo:

dove  e  rappresentano i pesi delle due variabili considerate.

Il valore dei pesi da assegnare si trova risolvendo il sistema qui di seguito:

Il cokriging è quindi utile quando la variabile principale è sottocampionata rispetto alle altre, come nel nostro caso in cui si ha che la variabile pozzi01” (principale), che rappresenta i valori di falda misurati al 2001, ha un terzo dei campioni della variabile pozzi91 (di appoggio).