Falda alluvionale nell’Alta valle del Tevere

Una tesi di laurea sulla ricostruzione tramite GIS della superficie piezometrica della falda alluvionale nell’Alta valle del Tevere

Analisi dati pozzi 2001

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

6. Analisi dei livelli di falda con approccio stocastico

6.4 Analisi dati pozzi 2001

Anche per l’analisi di questo dataset sono stati eseguiti gli stessi passi svolti per il dataset precedente. In questo caso però per l’analisi dei dati è stata usata sia la statistica univariata che quella multivariata in quanto la variabile del 2001 era sottocampionata rispetto a quella del Maggio ’91.

Per poter analizzare i dati sono stati per prima cosa, importati dentro R

> pozzi01 <getSites6sp(“pozzi_pri01”)

si sono quindi analizzati i dati grezzi;

> hist(pozzi01$LIV_PRIO1,nclass=20)

> plot(density(pozzi01$LIV_PRIO1))

figura625

figura626

La presenza sulla curva di densità di più massimi relativi potrebbe far pensare alla presenza di più popolazioni, in realtà questi massimi relativi sono dovuti alla scarsità di dati del dataset (solo 27 punti per tutta l’area coinvolta), che non riescono a coprire tutti gli intervalli di frequenza. Possiamo quindi affermare, essendo poco realistica la nascita di una nuova falda nell’arco di tempo di 10 anni, che il dataset si riferisce ad un’unica falda come per il dataset del Maggio ’91.

È stata poi verificata la presenza di eventuali outliers attraverso l’analisi del boxplot.

figura627

Come si può vedere dal grafico questo dataset presenta 4 outliers di cui soltanto uno è molto distante dal 95° percentile. È stato osservato che la presenza di questi outliers non è comunque legata ad errori grossolani di misura, si è deciso allora, anche in virtù della scarsa quantità di dati, di considerare anche questi punti per le operazioni di interpolazione.

6.4.1 Costruzione variogrammi

Si è quindi passati alla costruzione dei variogrammi:

> plot(variogram(LIV_PRIO1~1,loc=pozzi01,cutoff=8000,width=800))

figura628

La Figura 6.28 mette in luce, come per i dati precedenti, la presenza di un drift. Si è provveduto a detrendizzare il dato in modo da studiare solo la parte stocastica. Si sono quindi andati a studiare le profondità della falda dal piano campagna sostituendo ~1 con ~DEM.

> plot(variogram(LIV_PRIO1~DEM,loc=pozzi01,cutoff=8000

width=800))

figura629

Come si può vedere dalla Figura 6.29 il variogramma non risulta chiarissimo, presentando forti oscillazioni tra un lag e il successivo, questo perché i punti a disposizione sono pochi.

Sono poi state cercate eventuali anisotropie del dato, plottando la mappa del variogramma (Figura 6.30)

>plot(variogram(LIV_PRIO1~DEM,loc=pozzi01,cutoff=8000,

width=800,map=T))

figura630

Anche in questo caso si vede, anche se non chiaramente, come per i dati del ’91, che la direzione di maggior continuità è quella a 135° rispetto al nord.

Sono quindi stati analizzati i variogrammi per le direzioni 135° e 45°

>plot(variogram(LIV_PRIO1~DEM,loc=pozzi01,cutoff=8000,width=800

alpha=c(45,135)))

figura631

6.4.2 Modellazione variogramma

Come per i dati del ’91 i punti sono stati interpolati con due modelli esponenziali e due modelli sferici (manualmente e automaticamente) e calcolato il rapporto di anisotropia (dato dal rapporto tra i range per le due direzioni).

Inizialmente il variogramma è stato fittato in modo automatico.

>

plot(variogram(LIV_PRIO1~DEM,loc=pozzi01,cutoff=8000,width=800,alp

ha=c(135,45)),fit.variogram(variogram(LIV_PRIO1~DEM,loc=pozzi01,cu

toff=8000,width=800,alpha=c(135,45)),vgm(model=”Exp”,psill=36,nugg

et=5,range=3000,anis=c(135,0.6))))

>

plot(variogram(LIV_PRIO1~DEM,loc=pozzi01,cutoff=8000,width=800,alp

ha=c(135,45)),fit.variogram(variogram(LIV_PRIO1~DEM,loc=pozzi01,cu

toff=8000,width=800,alpha=c(135,45)),vgm(model=”Sph”,psill=36,nugg

et=5,range=3000,anis=c(135,0.6))))

Si è passato poi al fitting manuale

>

plot(variogram(LIV_PRIO1~DEM,loc=pozzi01,cutoff=8000,width=800,alp

ha=c(135,45)),vgm(model=”Exp”,psill=36,nugget=2,range=2000,anis=c(

135,0.6)))

>

plot(variogram(LIV_PRIO1~DEM,loc=pozzi01,cutoff=8000,width=800,alp

ha=c(135,45)),vgm(model=”Sph”,psill=30,nugget=2,range=3000,anis=c(

135,0.6)))

Per questo dataset il fitting manuale non è stato facile in quanto il numero ridotto dei punti campionati non ha permesso di avere dei variogrammi sperimentali molto chiari.

figura632

figura633

figura634

figura635

6.4.3 Cross-validation dei dati

Come già visto nel paragrafo 6.3.3 si esegue la cross-validation tramite il comando krige.cv, utilizzando come modelli di variogramma quelli appena descritti nella precedente sezione:

>

CV_fit_variogramma_pri01_exp=krige.cv(LIV_PRIO1~DEM,loc=pozzi01,mo

del=fit_variogramma_pri01_exp,nfold=nrow(pozzi01))

>

CV_variogramma_pri01_exp=krige.cv(LIV_PRIO1~DEM,loc=pozzi01,model=

variogramma_pri01_exp,nfold=nrow(pozzi01))

>

CV_fit_variogramma_pri01_sph=krige.cv(LIV_PRIO1~DEM,loc=pozzi01,mo

del=fit_variogramma_pri01_sph,nfold=nrow(pozzi01))

>

CV_variogramma_pri01_sph=krige.cv(LIV_PRIO1~DEM,loc=pozzi01,model=

variogramma_pri01_sph,nfold=nrow(pozzi01))

Analizziamo ora i residui della cross-validation per decidere il modello di variogramma migliore costruendo, gli istogrammi, la curva di densità e calcolando i principali indici statistici.

figura636

figura637

> summary(CV_fit_variogramma_pri01_exp$residual)

> summary(CV_variogramma_pri01_exp$residual)

> summary(CV_fit_variogramma_pri01_sph$residual)

> summary(CV_variogramma_pri01_sph$residual)

tabella3

L’analisi dei residui mostra come in questo caso si abbia un’accuratezza molto buona per quasi tutti i modelli (media prossima a zero) ma, una pessima precisione (RMSE abbastanza elevato) ed un range di valori molto elevato. Tutto questo porta ad affermare che la causa di questi risultati è la scarsità, in termini di quantità, del dataset.

È da notare inoltre come tutti i modelli tendano a sovrastimare i valori reali del dato in quanto presentano un MSDR>1.

6.4.4 Geostatistica multivariata sui dati del 2001

Al fine di ottenere un interpolazione più veritiera del dataset in questione (Primavera 2001), sono stati presi in considerazione anche i risultati ottenuti per l’altro dataset (Maggio ’91). Attraverso l’utilizzo della statistica multivariata si può interpolare un set povero di dati utilizzando oltre ai valori di questo anche quelli di un altro dataset.

Per lo studio del suddetto dataset oltre che al kriging universale si è quindi deciso di ricorrere al cokriging [D. G. Rossiter, 2006]. In questo modo nella fase di interpolazione si è tenuto conto del andamento del set di dati di appoggio (Maggio ’91) per la stima dei valori del dataset principale (Primavera 2001).

Per fare ciò tramite R bisogna prima di tutto creare un file che possa contenere tutti i dati necessari attraverso il comando gstat della omonima libreria. Il primo dato da inserire è la variabile da stimare (cioè i livelli della Primavera 2001).

> g=gstat(id=”pozzi01″,formula=LIV_PRIO1~DEM,location=pozzi01)

#crea il file oggetto g contenente i valori del file pozzi01, formula=LIV_PRIO1~DEM indica il variogramma sperimentale, id=”pozzi01″ assegna un identificativo all’interno di g; viene poi inserita la variabile di supporto dentro g.

>

g=gstat(g,id=”pozzi91″,formula=LEV_MAG91~DEM,location=pozzi91_meno

)

#aggiunge a g il variogramma sperimentale per il dataset pozzi91_meno;

Il file g, se chiamato, si presenta come in Figura 6.38

figura638

> plot(variogram(g,cutoff=8000,width=800))

#plotta i variogrammi per i due dataset contenuti in g ed un covariogramma calcolato sulla base dei due;

figura639

Come possiamo vedere il covariogramma segue lo stesso andamento dei due variogrammi che lo hanno generato, ciò significa che i due set di dati sono in qualche modo correlati tra loro.

Si è poi passati alla modellazione dei tre variogrammi andando ad inserire dentro g i tre modelli. L’unico vincolo matematico che si deve rispettare è quello che i tre variogrammi abbiano tutti lo stesso range. È stato di seguito riportato i comandi per il solo modello esponenziale perchè è stato osservato essere quello più aderente al covariogramma sperimentale.

>

g=gstat(g,id=”pozzi01″,model=vgm(model=”Exp”,nugget=5,psill=36,ran

ge=3000),fill.all=T)

#fill.all=T inserisce per i tre variogrammi lo stesso modello, quindi con lo stesso range;

Il passo successivo è il fitting dei tre modelli tramite il comando fit.lmc

> g=fit.lmc(variogram(g,cutoff=8000,width=800),g)

#fitta automaticamente i tre variogrammi utilizzando per tutti lo stesso range (lmc=linear model coregionalization);

> plot(variogram(g,cutoff=8000,width=800),model=g$model)

#plotta i tre variogrammi sperimentali e i tre modelli fittati (Figura 6.40);

figura640

Finita la modellazione dei variogrammi si è passati alla cross-validation dei dati.

Il comando che permette la cross-validation multivariata è gstat.cv

> CV_c_exp=gstat.cv(g,nfold=nrow(pozzi01))

#svolge la cross-validation della variabile pozzi01 togliendo un pozzo alla volta;

È stata successivamente svolta un’analisi dei residui ottenuti, plottando per prima cosa l’istogramma dei residui (Figura 6.41)

> hist(CV_c_exp$residual,nclass=20)

e facendo un summary e calcolando i valori di RMSE e MSDR (Tabella 4).

> summary(CV_c_exp$residual)

figura641

tabella4

Dai risultati ottenuti da questa cross-validation si vede chiaramente che il modello multivariato restituisce risultati migliori rispetto al modello univariato. Si ha infatti una maggiore precisione in quanto il valore RMSE è passato da ~7 a ~3,5.

L’accuratezza invece, data dalla media dei residui, è rimasta pressoché invariata. Il termine MSDR, che esprime il rapporto tra la varianza vera e quella calcolata nella crossvalidation, dice che i valori interpolati tendono ad essere sottostimati in quanto ha un valore minore dell’unità. Questo comportamento è un aspetto classico del kriging che tende a “lisciare” il comportamento di una variabile soprattutto se tale variabile risente molto di effetti di piccola scala.

Scelto quindi come modello per la costruzione della superficie della falda quello multivariato si procede alla costruzione della stessa.

> CK=predict.gstat(g,newdata=DEM)

Completato il cokriging la mappa può essere stampata a video attraverso il comando spplot:

> spplot(CK,zcol=1,main=”Altezze piezometriche Primavera 2001”)

ottenendo la mappa in Figura 6.42

Questa è stata infine importata in Grass (Figura 6.43) attraverso il comando:

> writeRast6sp(CK,”CK”,zcol=”pozzi01.pred”)

figura642

figura643

Una volta trovata la superficie piezometrica della primavera 2001 è stato possibile calcolarne la soggiacenza dal piano campagna.

> r.mapcalc “prof_falda_pri01=dtm25CK”

6.4.5 Analisi dati Luglio 1975 e Estate 1999

I dati relativi a queste due date (reperiti solo in corso d’opera) non sono stati trattati con la geostatistica (kriging) in quanto i risultati ottenuti con il RST non differiscono molto da quelli ottenuti con il kriging.

6.5 Confronti tra livelli di falda

Ricostruite le due superfici per i due dataset, l’ultimo passo è stato quello di confrontarle tra di loro, facendo una differenza tra le due superfici, per osservare eventuali variazioni del livello di falda (Figura 6.44).

> r.mapcalc “diff_UK_CK=UKCK”

figura644

La mappa diff_UK_CK delle differenze di quota piezometrica mostra come non ci siano grosse variazioni di quota della falda nel corso di questi 10 anni. L’unica zona in cui si hanno variazioni positive, quindi abbassamenti, di ~3/4 metri è nella zona nordovest subito a valle dell’invaso di Montedoglio. Ciò è forse dovuto alla fase di erosione e di approfondimento dell’alveo che subiscono i corsi d’acqua subito a valle degli sbarramenti, in quanto diminuisce la capacità di ricarica della falda da parte del fiume. Scendendo verso sudest si assiste ad una progressiva diminuzione delle differenze, fino ad arrivare nella zona di Città di Castello dove la mappa mostra addirittura innalzamenti di ~3/4 metri.

Autore: Penco100

Condividi questo articolo su

Invia commento

Il tuo indirizzo email non sarà pubblicato.