Interpolazione deterministica dei livelli di falda

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

5. Interpolazione deterministica dei livelli di falda

Un interpolatore deterministico si basa sull’interpolazione di dati puntuali tramite funzioni matematiche, come possono essere ad esempio le funzioni “spline”. Questi interpolatori assegnano un singolo valore per ogni punto della mappa, dato dalla media pesata dei valori circostanti noti. Esistono anche altri metodi di interpolazione deterministici, come:

• inverso della distanza pesata (IDW)
• regressioni polinomiali
• regressioni spline

5.1 Regularized spline with tension

Per la ricostruzione della falda tramite il metodo deterministico si è usato il modulo v.surf.rst di Grass, il quale utilizza come metodo di interpolazione la regularized spline with tension.

La regularized spline with tension è un metodo di interpolazione che approssima la superficie attraverso le funzioni spline, cioè funzioni costituite da un insieme di intervalli di polinomi di grado p raccordati tra loro. Lo scopo è interpolare in un intervallo circoscritto, un insieme di punti (detti nodi della spline), in modo da creare una funzione continua (almeno fino ad un dato ordine di derivate) in ogni punto dell’intervallo. Mentre l’interpolazione lineare utilizza una funzione appunto lineare per ciascuno degli intervalli, l’interpolazione spline si serve dei suddetti intervalli di polinomi di grado piccolo (massimo terzo grado), scegliendoli in modo che due polinomi successivi si saldino in modo liscio (ammette derivate di qualsiasi ordine). La funzione che si ottiene con un procedimento di questo genere si chiama funzione “spline”. La funzione interpolante ottenuta con la interpolazione spline (come quella ottenuta con la interpolazione polinomiale), rispetto alla interpolante ottenuta con l’interpolazione lineare, presenta errori inferiori ed è più liscia. L’interpolante spline risulta comunque più facile da valutare dei polinomi di grado elevato richiesti dall’interpolazione polinomiale, non soffre inoltre del fenomeno di Runge (problema relativo all’interpolazione polinomiale con polinomi di grado elevato), fenomeno per cui all’aumentare del grado del polinomio, l’interpolazione risultante oscilla in ampiezza verso gli estremi della funzione. Runge dimostrò che l’oscillazione poteva essere ridotta usando le funzioni spline, cioè suddividendo la curva in più parti in modo da utilizzare funzioni di grado inferiore.

Concettualmente, questo interpolatore assomiglia ad un foglio di gomma che passa attraverso i punti di input mentre minimizza la curvatura totale della superficie. Esso adatta una funzione matematica ad uno specificato numero di punti più vicini, mentre passa attraverso i punti di campionamento. Questo metodo è un metodo di interpolazione locale, cioè opera nell’intorno del punto da interpolare ed è abbastanza buono quando si devono interpolare campioni che variano dolcemente nello spazio. Un difetto di questo interpolatore (deterministico) è però l’impossibilità di quantificare l’errore commesso.

I parametri su cui si può lavorare per ottenere un approssimazione migliore sono la tension e lo smoothing. Per determinare i parametri che produrranno una superficie con le proprietà volute, è utile sapere che il metodo è dipendente dalla scala e la tension funziona come parametro di rescaling, cioè tension alte aumentano le distanze tra i punti e riducono il range d’influenza di ogni punto; tension basse diminuiscono le distanze tra i punti e aumentano il range d’influenza. Il parametro smoothing regola lo scostamento tra il valore reale del dato e quello calcolato, più specificatamente per un valore di smoothing uguale a zero si ha che la superficie interpolata passa esattamente per il valore reale in quel punto; al crescere di questo valore gli scostamenti tra il valore vero e quello interpolato aumentano producendo un “lisciamento” della superficie. La tension ha un importante effetto sulla superficie interpolata: la fa comportare come una sottile membrana, quindi producendo picchi e crateri, in corrispondenza dei punti campionati, per valori di tension elevati; oppure come una lastra di acciaio (poco flessibile), per bassi valori di tension. I parametri ottimali di tension e smoothing possono essere determinati tramite la procedura di crossvalidation (vedi paragrafo 5.2) dei dati; per riuscire a determinare la giusta combinazione si deve eseguire la procedura in modo iterativo, modificando un parametro alla volta fintanto che non si minimizza l’errore statistico sul residuo (valore verovalore stimato).

La funzione è data dalla somma di una “trend function” ed una “radial basis function”:

dove z (r) il valore interpolato,  è il peso j – esimo e  è funzione di base radiale con una particolare forma che dipendente dalla scelta dei pesi.

La trend function è data da:

dove f l (r) è un sistema di variabili indipendenti.

5.2 Cross-validation

La cross-validation è un modo per verificare le assunzioni del modello usato. Il metodo di funzionamento è molto semplice, consiste nell’estrarre ogni campione a turno, assumendolo come dato mancante, per poi fare la stima, nello stesso punto, del valore in base ai dati che gli stanno attorno. Quindi per ogni campione abbiamo:

• il suo valore vero
• il suo valore stimato

Dopo di che i valori stimati e i valori reali vengono paragonati. La loro differenza è detto “residuo della cross-validation”. Il residuo sarà poi studiato statisticamente per verificare l’attendibilità del modello scelto. Dai dati in uscita dalla cross-validation si possono calcolare gli errori sperimentali e gli errori previsti, attraverso le quantità medie: errore medio e varianza dell’errore medio; questi due valori possono anche essere standardizzati, ottenendo l’errore medio standardizzato e la varianza dell’errore medio standardizzato. I valori di questi parametri danno il grado di correttezza del modello. La cross-validation ci permette quindi di paragonare l’impatto dei differenti modelli sui risultati dell’interpolazione.

5.3 Applicazione ai dati

Dopo aver preparato i set di dati (Luglio 1975, Maggio 1991 Giugno 1999 e Primavera 2001), gli stessi sono stati elaborati dal modulo v.surf.rst (vedi sezione 5.4), per ricostruire la superficie piezometrica della falda in questione per i suddetti periodi. Questo modulo è utilizzato soprattutto per la costruzione dei DEM in quanto è un ottimo interpolatore per questo tipo di superfici.

v.surf.rst permette di eseguire l’interpolazione di punti disposti nello spazio basandosi sulle loro coordinate z o sui loro attributi tabellari (in questo caso gli attributi sono le quote rispetto al livello del mare della falda). I dati in ingresso da inserire per l’elaborazione sono i vettoriali dei pozzi già descritti, in uscita il modulo offre diversi output in quanto permette la creazione di diverse mappe raster, tra cui il DEM, la carta delle pendenze e quella delle superfici; si ha anche la possibilità, prima di creare queste mappe, di fare una cross-validation dei dati e salvarla in un file vettoriale, che poi può essere analizzato statisticamente, tramite il modulo v.univar, per verificare l’attendibilità dei risultati ottenuti. v.univar restituisce come output il valore medio dell’errore, la deviazione standard ed il range di valori, che permettono di capire quale sia la migliore combinazione dei parametri tension e smoothing.

Sono state eseguite le cross-validation per le superfici interpolate relative ai diversi dati a disposizione e si è verificato che i parametri per cui si hanno i risultati migliori, in termini di media e varianza dell’errore, sono stati:

1. tension uguale a 80 e uno smoothing uguale a zero senza nessuna anisotropia per la superficie del Luglio 1975
2. tension uguale a 40 e uno smoothing uguale a zero senza nessuna anisotropia per la superficie del Maggio 1991;
3. tension uguale a 40 e uno smoothing uguale a zero senza nessuna anisotropia per la superficie del Giugno 1999
4. tension uguale a 60 e una smoothing uguale a zero senza nessuna anisotropia per la superficie della Primavera 2001.

Gli errori commessi (ottenuti tramite cross-validation) utilizzando questi parametri sono riassunti in Tabella 1.

Le superfici finali ottenute, per il 1991 e il 2001, sono rappresentate in Figura 5.1 e 5.2

Si osserva chiaramente che le quote della falda mostrano trend che va da NW verso SE seguendo l’andamento delle quote del territorio (come era intuibile)

L’interpolazione per la Primavera ’01, è risultata molto più approssimativa a causa della minore quantità di punti (solo 27 contro i 87 del Maggio ’91).

Si sono quindi confrontate le due superfici facendo una semplice differenza delle altezze piezometriche per le due date tramite il modulo r.mapcalc, il quale permette di eseguire comandi di “map algebra” sui raster (addizione, sottrazione…).

Interrogando la mappa delle differenze piezometriche (Figura 5.3), tramite l’apposito comando di query, si è visto che:

1. nei punti in cui si hanno pozzi campionati nel 1991 e nel 2001 le differenze piezometriche non superano quasi mai il metro;
2. nei punti in cui non si hanno pozzi comuni si ottengono invece forti differenze.

Queste differenze di livello, non devono trarre in inganno. È infatti impensabile avere punti in cui la falda si è abbassata di 60 metri (colore blu e rosso) e punti in cui si è alzata di 15 metri (zona gialla). Riguardo alla zona di ipotetica depressione, queste grosse differenze sono dovute sicuramente al fatto che lì si hanno pozzi campionati solo per l’anno ’91, mentre nel ’01 il pozzo più vicino a questa zona è a quasi 2 Km di distanza ed a 50 metri di dislivello.

Nelle zone in cui si ottengono presunti innalzamenti di 15 metri (zone gialle) possiamo ipotizzare che questi siano legati al fatto che, durante le misure del ’91, si stavano eseguendo degli emungimenti che potrebbero aver provocato una temporanea depressione della zona. Oppure, visto che nelle zone gialle non troviamo mai un pozzo del ’01, potrebbe essere, come per la zona precedentemente analizzata, un errore di interpolazione.

Per avere un idea ancora più completa della evoluzione della falda, tramite il modulo v.surf.rst, è stato ricostruito il livello di falda per l’estate 1975. La superficie ottenuta dall’interpolazione dei punti (218 pozzi) è stata quindi confrontata con la superficie del Estate 1999 (28 pozzi). Anche questo confronto non ha portato a differenze significative tra le due superfici piezometriche come mostra la Figura 5.4.

Le differenze più forti (maggiori di 3 metri) tra le due superfici si trovano solo ai bordi e nelle zone non coperte da pozzi in entrambe le date (zona gialla, blu e rossa). Ciò significa che queste differenze sono dovute ad errori di interpolazione della superficie e non ad effettive depressioni della falda.

L’interpolatore spline segue, come detto precedentemente, un modello deterministico di interpolazione quindi non tiene in considerazione l’incertezza associata alle variabili di input (al contrario dei modelli stocastici come il kriging). È importante sottolineare che ciò che differenzia i due metodi è che in quello stocastico, si tiene conto della variabilità dei dati in input. Generalmente questo metodo ha una struttura più complessa di quello deterministico e presenta una maggiore complessità nei calcoli. I modelli stocastici sono solitamente molto più affidabili e forniscono risultati più aderenti alla realtà.

Per questo motivo si è cercato, come descritto nel prossimo capitolo, di studiare il fenomeno con un metodo stocastico come il kriging.

5.4 Procedura creazione superfici con v.surf.rst

Si riportano qui di seguito alcuni passaggi relativi alla procedura di interpolazione descritta nella sezione precedente.

Partendo dai vettoriali contenenti i dati si sono digitate le seguenti stringhe di comando (su Grass GIS):

-v.surf.rst inpunt=pozzi91 layer=1 zcolumn=LEV_MAG91

cvdev=cv_t40_s0_noan maskmap=MASK zmult=1,0 Tension=40 smooth=0 -c

#fa la crossvalidation dei dati per i valori di tension=40, smoothing=0 e senza anisotropia (input=vettoriale d’ingresso, zcolumn=colonna con i valori da interpolare, cvdev=file di uscita

della cross-validation dei punti, maskmap=raster usato come maschera, zmult=fattore di conversione dei dati, tension=valore di tension, smooth=valore dismoothing, -c=bandierina per fare la cross-validation).

– v.univar map=cv_t40_s0_noan type=point column=flt1 layer=1

#fa l’analisi statistica dei residui salvati nel file di input (map=file di ingresso da analizzare, type=tipo di input, column=colonna dove risiedono i dati)

queste due stringhe sono state ripetute per vari valori dei parametri tension e smoothing fino ad ottenere gli errori più bassi. Trovati i valori che minimizzano gli errori si è passati alla creazione della superficie interpolata.

– v.surf.rst inpunt=pozzi91 layer=1 zcolumn=LEV_MAG91

elev=LEV_MAG91 maskmap=MASK zmult=1,0 Tension=40 smooth=0

#questo comando produce in uscita una mappa raster della superficie del livello della falda del Maggio ’91 (Figura 5.1)

– v.surf.rst inpunt=pozzi01 layer=1 zcolumn=LIV_PRIO1

elev=LEV_PRI01 maskmap=MASK zmult=1,0 Tension=60 smooth=0

#questo comando produce in uscita una mappa raster della superficie del livello della falda della Primavera ’01 (Figura 5.2)

Una volta trovate le due superfici piezometriche si è prodotta la mappa delle differenze piezometriche

– r.mapcalc “differenze_91_01=LEV_MAG91LEV_PRI01”

#questo comando produce una mappa contenente in ogni singolo pixel le differenze tra il livello del ’91 e quello del ’01 dello stesso pixel (Figura 5.3)

Si è poi creato il vettoriale DIFF91_01 delle curve di livello del raster differenze_91_01 per avere una visione più chiara della superficie

– r.contour input=differenze91_01 output=differenze91_01

minlevel=-20 maxlevel=75 step=5

#questo comando ha prodotto un vettoriale contenente le curve di livello a partire dal raster differenze91_01 (input=raster in input, output=vettoriale di uscita, minlevel=valore minimo della isolivello, maxlevel=valore massimo della isolivello, step=passo isolivello).

Le operazioni sopra citate sono state eseguite anche per i dati del 1975 e del 1999.

Tutte le operazioni svolte da riga di comando possono anche essere eseguite tramite una semplice interfaccia grafica (Figura 5.5) che permette di scrivere i comandi semplicemente compilando i campi necessari.