Università degli studi di Padova - Metodologie geostatistiche...

1

Sede Amministrativa: Università degli Studi di Padova

Dipartimento di Geoscienze

Dottorato di ricerca in Scienze della Terra

Ciclo XXIII

Metodologie geostatistiche nel idrogeologico

Direttore della Scuola: Ch.mo.Prof. Gilberto Artioli

Supervisore : Ch.mo Dr. Paolo Fabbri

Dottoranda : Dr. Amany Hammam

2

INDICE

Pagina

ABSTRACT. 4

SOMMARIO 6

CAPITOLO 1

1.1. INTRODUZIONE 7

CAPITOLO 2

2.1. INTRODUZIONE 11

2.2. GEOSTATISTICA 12

2.3. LE VARIABILI REGIONALIZZATE (VR) 13

2.4. VARIABILI ALEATORIE 14

2.5. MOMENTI 18

2.6. CONTINUITA SPAZIALE 22

2.7. VARIOGRAMMA 30

2.8. INTERPOLAZIONE SPAZIALE 35

2.9. METODI DETERMINISTICI 36

2.10. METODI STOCASTICI 38

2.11. KRIGING 39

2.11.1. KRIGING SEMPLICE 41

2.11.2. KRIGING ORDINARIO 47

2.12. COKRIGING 49

3

2.13. KRIGING INDICATORE 55

2.14. CROSS VALIDATION 58

CAPITOLO 3. Il software R


3.2. CHE COS‟è R? 61

3.3. LE FUNZIONALITÀ DELL‟AMBIENTE R 62

3.4. ALCUNI PACKAGES DI R 63

3.5. PACKAGE GSTAT 63

CAPITOLO 4.


4.2. LE ACQUE SOTTERRANEE 68

4.3. LA PIANURA VENETO 70

4.3.1. L‟ALTA PIANURA 71

4.3.2. MEDIA PIANURA 71

4.3.3. BASSA PIANURA 71

4.4. ANALISI GEOSTATISTICA 74

4.5. ANALISI DEI RISULTAT 81

4.6. CROSS VALIDAZIONE DEL KRIGING 91

4

CAPITOLO 5.


5.2. LE ACQUE SOTTERRANEE 97

5.3. RAPPRESENTAZIONE SPAZIALE DEI DATI 98

5.4. ANALISI GEOSTATISTICA 100

5.5. ANALISI DEI RISULTATI 104

5.5.1. KRIGING ORDINARIO 104

5.5.2. COKRIGING 109

4.6. CROSS VALIDAZIONE 107

CONCLUSIONE 111

BIBLIOGRAFIA 113

5

ABSTRACT

This thesis is concerned with environmental and hydrogeological issues

approached by geostatistical methods. The objectives of this study are to apply

geostatistical tools (like variogram, kriging methods) to different hydrogeological

parameters for mapping the data.

The first part of the work is dedicated to presenting the geostatistical

theory, giving more weight to those aspects (exploratory data analysis, indicator

coding of information, using of secondary information) that play an important role

in hydrogeological and environmental issues.

Then, two hydrogeological and environmental issues, whose sites are all

located in Veneto Region (Italy), are approached by using geostatistical methods.

Case Study 1. Spatial distribution of Ground Water Temperature in Veneto

Region

The geostatistical methodology is applied to a groundwater Temperature

case in an area situated near Padova (northeast Italy). The analysis of 145

ground water temperature measurements, ranging between 10°C and 17°C, and

of 100 meters in depth. This pattern in the spatial distribution of temperatures to

some extent seems related to the structural framework of the area and should be

taken into account when quantifying spatial distribution of temperature.

Geoatatistical tools including Ordinary Kriging and Indicator kriging are applied

for ground water temperature values.

6

Case study 2. Spatial distribution of Chlorine in the Geothermal Euganean

field

The Euganean geothermal reservoir is principally located in fractured

Mesozoic limestone. The analysis of 19 ground water Chlorine measurements,

200 meters in depth, ranging between 1410 mg/l and 2750 mg/l, reveals that the

spatial distribution of Chlorine along the geothermal area is characterized by

some degree of spatial continuity. Ground water Temperature measurement

ranged between 60°C and 86°C. Finally, results are presented and discussed for

Ordinary Kriging and Cokriging.

7

SOMMARIO

In questa tesi vengono affrontate una serie di problematiche di tipo

idrogeologico ed ambientale secondo la metodologia geostatistica. Gli obiettivi di

questo studio sono quelli di applicare gli strumenti della geostatistica

(variogrammi, metodi kriging, etc.) a deversi parameteri idrogelogici per la

mappatura dei dati.

Nella parte iniziale della tesi, ci si dedica all‟esposizione e

all‟approfondimento di queste tematiche (analisi esplorativa dei dati, ecodifica

indicativa dell‟informazione) che risultano particolarmente rilevanti nel contesto

idrogeologico ed ambientale.

Caso Studio 1. Distribuzione spaziale delle temperatura delle acque

sotterranee nella Regione Veneto

I metodi di analisi geostatistici sono stati applicati allo studio delle

temperatura delle acque sotterranee, nella Region Veneto (Italia). L‟analisi dei

145 valori di temperatura delle acque sotterranee, oscillanti tra i 10 ºC e gli 17 ºC,

ad un profondità di 100 m. La struttura di continuità spaziale sembra in qualche

modo relazionata alle caratteristiche tettonico-strutturali dell‟area, questo fattore

deve essere tenuto in considerazione nel processo di distribuzione spaziale delle

temperature. Sono stati applicati per i valori della temperatura delle acque

sotterranee il Kriging Ordinario ed il Kriging indicatore.

Caso studio 2. Distribuzione del Cloruri nel bacino Termale Euganeo

Il reservoir termale considerato e` localizzato principalmente nel calcare

fratturato mesozoico dell‟area euganea (Regione Veneto, Italia). L‟analisi di 19

valori di Cloruri delle acque termali, oscillanti tra i 1410 mg/l e gli 2750 mg/l, ad

un profondità di 200 m, evidenziano la presenza di una certa continuità nella

distribuzione spaziale dei Cloruri. Vengono inoltre confrontati i valori di cloruri con

i valori di temperatura delle acque termali, oscillanti tra i 60 ºC e gli 86 ºC.

Conseguentemente, è stato sffettuata l‟analisi attraverso il Kriging Ordinario e

Cokriging.

8

CAPITOLO 1

1.1. INTRODUZIONE

Questa tesi è nata con lo scopo di fornire una metodologia per la

modellazione e la previsione spaziale della concentrazione del cloruri e delle

temperature delle acque sottorranee, a partire da misure effettuate in alcuni pozzi

campionati. Attraverso l‟utilizzo di un potente software, quale è il linguaggio R

(Versione 2.6.2), è stato possible condurre l‟intera analisi, dall‟esplorazione

preliminare alla previsione spaziale di specifici elementi a determinate profondità.

Con altri software è, invece necessario spezzare le fasi di studio, importando i

dati da un ambiente all‟altro, e quindi vengono enormemente dilatati i tempi di

raggiungimento dei risultati.

Lo sviluppo della geostatistica è iniziato negli anni ´50 del secolo scorso

quando, in Sud Africa, D. Krige cominciò ad applicare le tecniche statistiche per

le previsioni delle riserve minerarie. A partire degli anni ´60 il matematico

francese G. Matheron formulò le basi teoriche di tali metodi e un modello che

tenesse conto delle fluttuazioni delle variabili osservate ed, in particolar modo,

della correlazione spaziale.

Matheron (1971) ha definito le variabili spaziali come variabili

regionalizzate per descrivere la loro natura, espressa da due aspetti

complementari: uno casuale, che tiene conto delle irregolarità locali presentati nel

fenomeno (piccola scala), e un aspetto strumentale, che riflette le caratteristiche

globali dello stesso (larga scala).

9

Le variabili regionalizzate sono la realizzazione di processi stocastici, e

uno degli obiettivi della geostatistica è quello di identificare le due componenti

che le caratterizzano, a partire dalle informazioni fornire dal campione osservato.

Nella maggiore parte dei casi è impossible controllare l‟assunto che una variabile

sia la realizzazione di una funzione casuale, avendo a che fare con una singola

realizzazione.

E` necessario allora imporre delle ipotesi di stazionarietà, tramite le quali è

possible ritenere il fenomeno omogeneo, e si può definire la funzione

variogramma che in geostatistica descrive la correlazione spaziale. Dagli anni ´50

la geostatistica è un mezzo per descrivere la variabilità spaziale e l‟applicazione

più diffusa, il Kriging, consiste nella previsione dei valori di una grandezza in

punti in cui non è possible effettuare misurazini. Tale obiettivo si può raggiungere

utilizzando le osservazioni disponibili nei punti campionati, supponendo

l‟esistenza di una relazione tra valori osservati e quelli incogniti. Proprio per

questo la Geostatistica, come detto sopra, è stata inizialmente utilizzata

dall‟industria minieraria per ridurre gli elevati costi delle perforazioni. In seguito si

sono rivelati i molteplici campi di applicazione di tale disciplina: scienze del suolo,

idrogeologia, meteorologia, scienze ambientali, agricoltura, ingegneria civile, ecc.

L‟applicazione della geostatistica ai problemi ambientali (ad esempio, le

acque sotterranee contaminante) si è anche dimostrato un potente strumento di

previsione, consentendo il coordinamento delle attività di acquisizione dei dati per

progettare sistemi di monitoroggio. L‟a raccolta dei dati è spesso incompleta, con

conseguente incertezza nella comprensione del problema. Questa incertezza

crea il problema del campionamento supplementare, e di conseguenza

aumentano di molto i costi di analisi . La Geostatistica offre un mezzo per

quantificare questa incertezza, sfruttando i dati esistenti a sostegno di

campionamento ottimale.

10

Le procedure utilizzate per l‟analisi geostatistica coinvolgono sia tecniche

puramente statistiche, sia procedimenti cartografici per mezzo dei quali si ha un

riscontro visivo dei risultati matematici ottenuti e si può comprendere in modo

migliore i comportamenti del fenomeno studiato. Spesso si hanno a disposizione

quantità limitate di osservazioni, riguardo ad una variabile di interesse, con le

quali costruire le mappe di previsione, e perciò queste avranno un basso grado di

risoluzione e un‟elevata varianza di previsione. In situazioni di questo tipo si

ricorre a informazioni supplementari quali possono essere as esempio la

stratigrafia o la conoscenza delle attività umane svolte nella zona di interesse, nel

caso si eseguano studi di idrogeologia. Gli obiettivi di questo studio sono di

applicare la geostatistica e sturmenti come variogrammi e metodi kriging) a

deversi parameteri idrogelogici per la mappatura dei dati.

11

CAPITOLO 2

2.1. INTRODUZIONE

In ambito geologico, oltre ad avere una grande importanza il valore

numerico del parametro da studiare, è importante la posizione che il dato

possiede nello spazio. I metodi statistici classici non tengono conto

dell‟informazione spaziale che invece è propria del dato geologico. La

geostatistica quindi offre un modo per descrivere e sfruttare l‟informazione

spaziale, permettendo inoltre di studiare la continuità spaziale che è un aspetto

essenziale di molti fenomeni geologici e naturali in genere.

La geostatistica studia i fenomeni naturali che si sviluppano su base

spaziale a partire da informazioni derivanti da un loro campionamento. In

particolare, studia la variabilità spaziale dei parametri che descrivono i suddetti

fenomeni estraendone le regole in un quadro modellistico di riferimento e

usandole per effettuare le operazioni volte a dare soluzione a specifiche

problematiche riguardanti la caratterizzazione e la stima dei fenomeni stessi.

In tutti quei settori delle scienze in cui i fenomeni di studio hanno carattere

spaziale. In relazione alle applicazioni registrate negli ultimi tre decenni, tra i

settori applicativi si possono citare: le scienze geologiche e minierarie, l‟idrologia,

l‟idrogeologia, la scienza dei suoli, l‟agronomia, la geotecnica, la geofisica, il

telerilevamento, la climatologia, la meterologia, l‟oceanografia, le scienze

forestali, la zoologia, l‟epidemiologia, l‟igiene ambientale, ecc. Queste sono tutte

discipline che si collocano in un contesto di analisi dei dati di tipo spaziale o

spaziotemporale, e che sono alla base dello studio dei fenomeni ambientali che

più frequentemente si manifestano sul territorio (Raspa, 2000).

12

2.2. GEOSTATISTICA

La tipologia dei fenomeni ambientali è molto ampia. Essi vanno da quelli

più semplici, che si manifestano nello spazio e sono caratterizzati da una sola

variabile, a quelli più complessi, regolati da un gran numero di variabili, che

interagiscono tra loro e che sono definite non solo nello spazio, ma anche nel

tempo.

La geostatistica è una branca della statistica applicata che studia i

fenomeni che si manifestano su base spaziale. In particolare studia la variabilità

spaziale dei parametri che caratterizzano il fenomeno stesso, estrapolando un

modello matematico di riferimento (Raspa, 1995). Con “geostatistica” si intende

un insieme di metodi statistici e probabilistici che sono stati sviluppati

dall‟industria mineraria in Sud Africa negli anni cinquanta (Krige, 1951) a partire

da analisi fatte in Scandinavia nel settore forestale. I principi teorici di questa

recente disciplina sono stati sviluppati da Matheron (1963), chiamandoli “teoria

delle variabili regionalizzate”.

Le operazioni di tipo geostatistico sono fondamentalmente dirette alla

caratterizzazione spaziale o temporale di fenomeni fisici e alla loro

modellizzazione e stima risultati delle operazioni, oltre a rappresentare elementi

per l‟interpretazione e comprensione dei fenomeni, trovano impiego nella messa

a punto dei sistemi di supporto alle decisioni e nella progettazione e gestione dei

sistemi di monitoraggio.

In molte applicazioni che riguardano l‟igiene ambientale, spesso si pone il

problema di stimare la probabilità che la concentrazione di un inquinante, o di un

parametro da tenere sotto controllo, superi, in un punto del tempo e dello spazio,

un determinato valore di soglia.La geostatistica è preferita alle tecniche della

statistica classica per due ragioni: 1. la statistica classica considera le

osservazioni indipendenti e non tiene conto della distribuzione spaziale tipica

della maggior parte dei fenomeni ambientali (Oliver et al., 1989); 2. la

geostatistica utilizza le informazioni spaziali contenute nei dati per stimare il

valore della variabile in aree non campionate (Issaks e Srivastava, 1989).

13

La parte teorica di seguito e la simbologia relativa alle formule è stata

tratta da “An introduction to applied Geostatistics” (Isaaks e Srivastava, 1989) e

“Geostatistical Software Library and User‟s Guide” (Deutsch e Journel,1998).

In geostatistica si introduce il concetto di “fenomeni regionalizzati” per

indicare quei fenomeni naturali che si concretizzano nello spazio e nel tempo

manifestando un determinato comportamento. Le variabili che determinato questi

fenomeni vengono individuale come delle variabili a scala regionale o variabile

regionalizzate.

2.3. LE VARIABILI REGIONALIZZATE (VR)

Un fenomeno naturale distribuito nello spazio è caratterizzato da una

distribuzione spaziale di un certo numero di quantità misurabili chiamate Variabili

Regionalizzate. (ad esampio il tenore di un minerale nel sottosuolo, la

temperatura dell‟aria in una stanza, ecc.). Una variabile regionalizzata è la

realizzazione di un processo stocastico, cioè di un insieme di variabili casuali

corrispondenti ai punti di un dominio 𝐷 oggetto di studio. Dunque, per ogni punto

𝑧(𝑥), in uno spazio 𝑑 dimentionale, il valore del parametro di interesse 𝑧(𝑥) è

una delle realizzazioni del processo 𝑍(𝑥).

Matheron (1971) ha definito le variabili spaziali come variabili

regionalizzate per descriverne la loro natura, espressa da due aspetti

complementari: uno casuale, che tiene conto di irregolarità locali presentate dal

fenomeno (piccola scala), e un aspetto strutturale, che riflette le caratteristiche

globali dello stesso (larga scala). Il modello fondamentale per una variabile

regionalizzata 𝑍(𝑢), osservata nella localizzazione 𝑥𝑖 , è:

𝑍 𝑥 = 𝑚 𝑥 + 𝑌 𝑥

La componente casuale è descritta da 𝑌 𝑥 e quella strutturale da 𝑚 𝑥 ,

chiamata trend. Nel caso in cui la componente casuale fosse assente, il

fenomeno potrebbe essere descritto dal solo trend, rappresentandolo tramite una

funzione.

14

Uno degli obiettivi della Geostatistica è quello di identificare le due

componenti suddette per risolvere il problema della previsione di una variabile

regionalizzata, a partire dall‟informazione frammentaria contenuta nel campione

osservato. L‟interpretazione dei parametri supporta l‟ipotesi che non sia possibile

descrivere completamente i fenomeni naturali attraverso regole deterministiche.

Nella maggior parte dei casi è impossibile controllare l‟assunto che il

parametro sia la realizzazione di una funzione casuale (come richiesto dalla

statistica classica) avendo a che fare con una singola realizzazione dalla

variabile.

2.4. VARIABILI ALEATORIE

Uno dei concetti fondamentali della teoria della probabilità è quello di

variabile aleatoria; prima di darne una definizione formale illustriamo questo

concetto con esempi.

Il numero delle particelle cosmiche che investono qualche porzione

della superficie terrestre in un dato intervallo di tempo è soggetto a

variazioni apprezzabili a causa di una moltitudine di fattori casuali.

La velocità di una molecola di gas varia nel tempo in conseguenza, ad

esempio, a collisioni casuali con altre molecole.

Questi esempi di numeri aleatori, che intervengono nei più diversi campi della

scienza e della tecnologia, conducono alla considerazione di quelle che vengono

chiamate “variabili aleatorie”. Osserviamo che, nonostante la diversità del

contenuto reale degli esempi precedenti, essi presentano nella sostanza

caratteristiche comuni dal punto di vista matematico.Infatti, in ciascun esempio

interviene una quantità che in un modo o nell‟altro descrive i fenomeni che si

stanno studiando. Sotto l‟effetto di circostanze casuali, ciascuna di queste

quantità può assumere una varietà di valori diversi e non si può stabilire in

anticipo quale valore la quantità assumerà, perche essa varia in modo casuale

da “prova” a “prova”. Solo dopo aver fatto la “prova” si potrà dire quale valore ha

assunto la variabile aleatoria.

15

Quindi, al fine di esaminare una variabile aleatoria, è necessario

conoscere tutti i valori (determinazioni) che essa può assumere. Ma una

semplice lista dei valori della variabile non è sufficiente per trarre delle

conclusioni sul fenomeno in esame. Infatti per specificare che cosa sia una

variabile aleatoria è necessario conoscere non solo il valore che essa può

assumere, ma anche quanto spesso, cioè con che probabilità essa assume

questo valore. La diversità delle variabili aleatorie è estremamente grande.

Il numero dei valori reali assunti da una variabile aleatoria può essere

definito come finito, più che numerabile. Introducendo dapprima la definizione di

variabile aleatoria di tipo finito per poi passare a una definizione più generica di

variabile aleatoria:-

Una variabile aleatoria x di tipo “finito” è una quantità variabile che può

assumere valori reali x1, x2, ......xn a seconda che, in dipendenza del caso, si

presenti uno degli eventi E1, E2, ......En , di probabilità note p1, p2, ......pn,

soddisfacenti la condizione: p1+p2, +......pn =1.

Dunque la variabile aleatoria di tipo “finito” nella forma canonica, è definita

come una quantità variabile che, in dipendenza del caso, può assumere valori

reali, tutto distinti, x1, x2,.....xn (determinazioni) a seconda che si presenti uno

degli incompatibili con probabilità rispettive p1, p2, ......pn, soddisfacenti la

condizione:

X1, X2 𝑝𝑖 = 1

𝑛

𝑖=1

Osserviamo che, in base a tale definizione, gli eventi suddetti E1, E2,

......En , costituiscono una base di eventi per l‟esperimento dalla suddetta variabile

aleatoria.

In generale si indica la variabile aleatoria con lettere maiuscole del tipo X,

Y, Z.......e le sue possibili determinazioni con le corrispondenti lettere minuscole

x1,..... xn o y1,....... yn o z1, ........z2 o ...... .

Vediamo ora un esempo di variabili aleatorie del tipo finito:-

16

Supponiamo di fare 100 misure di una data quantità (per es. la

temperatura di una zona del Nord Europa) e di ottenere dei valori cosi riportati:-

11 volte si ha la misura -2

15 volte si ha la misura -1

25 volte si ha la misura 0




Costruiamo la seguente tabella: Nella 1a riga compaionoi valori che la variabile

temperatura può assumere, nella 2 a riga la probabilità (desunta a posteriori) che

tali valori vengono assunti della variabile aleatoria.

La temperatura rappresenta in questo caso unaquantità che può assumere

vari valori, in dipendenza da circostanze metereologiche ritenute casuali. Di tale

variabile aleatoria infatti noi non conosciamo a priori un valore bene determinto,

ma tutti i valori numerici casuali che essa può assumere in

quell‟esperimento.L‟informazione però è incompleta se non specifichiamo anche

la probabilità che la variabile aleatoria assuma uno di quei valori. Nella 2a riga

scriviamo quindi tali probabilità: in tal modo la tabella rappresenta la “legge di

distribuzione” di quella variabile aleatoria.

In generale: siano xi (i = 1, ....... n) i valori tutti distinti che la variabile

aleatoria può assumere e p(xi) = pi, le rispettive probabilità (valutate a priori o a

seconda del problema).

si può affermare che la tabella:

x1 x2 xn

p1 p2 pn

rappresenta “la legge di distribuzione di una variabile aleatoria” se:

a) xi , xj P(xi e xj ) = 0

b) Almeno uno degli xi (i = 1.............n) è assunto della variabile aleatoria cioè

P(x1 o x2 ......... o xn ) = 1.

Le proprietà a) e b) possono essere sintetizzate nella relazione seguente:

17

𝑃𝑖

𝑛

𝑖=1

= 1

Una variabile aleatoria v.a. (o variabile casuale v.c. o variabile stocastica)

X è dunque una funzione misurabile a valori reali definita sullo spazio campione

Ω. Si definisce una v.a. quando si crea una corrispondenza fra insieme dei

risultati di una prova e insieme dei numeri reali (corrispondenza non

necessariamente biunivoca). In generale, possiamo dire anche che la variabile

regionalizzata è detta essere una realizzazione delle funzione aleatoria (figure

2.1 e 2.2).

Fig. 2.1: Variabile Regionalizzata (VR) come realizzazione della Funzione Aleatoria (FA)

Fig. 2.2: Il valore misurato in ogni punto è uno dei possibili risultati di una realizzazione della funzione aleatoria

18

Una variabile aleatoria si dice discreta se può assumere un numero finito o

un‟infinità numerabile di valori; si dice, invece, continua se può assumere i valori

in tutto R o in un suo sottoinsieme continuo. Una variabile aleatoria discreta

corrisponde ad un carattere quantitativo discreto, una variabile aleatoria continua

ad un carattere continuo. per una variabile aleatoria X, in analogia al caso

statistico, si definisce la funzione di ripartizione che ad ogni valore x, assumibile

da X, associa la probabilità che X sia non superiore a x: .

I principali parametri di una V.A. sono:

1- Sia X una variabile aleatoria discreta con un numero finito di valori;

ai valori x1, x2, x3, ......, xn corrispondono rispettivamente le

probabilità p1, p2, p3, ......, pn . Si dice valore atteso, o speranza

matematica, o semplicemente media di X la media aritmetica

ponderata di tali valori (assumendo le probabilità come pesi),

ovvero il valore E(X):

𝑬 𝑿 = 𝒙𝒊 𝒑𝒊

𝒏

𝒊=𝟏

La media di una variabile aleatoria, se esiste, gode di alcune

proprietà:

a) La media di una variabile aleatoria costante, è la costante

stessa.

b) Se X è una variabile aleatoria e a e b sono delle costanti allora:

E(aX+b) = aE(X)+b.

c) La media della somma di due variabili aleatorie è uguale alla

somma delle medie delle singole variabili. Cioè: E(X + Y) = E(X)

+ E(Y).

d) Sia data una variabile aleatoria X, consideriamo la differenza tra

la variabile aleatoria X e la sua media: X - E(X).

Questa nuova variabile aleatoria si chiama deviazione, o

variabile aleatoria scarto.

Si dimostra subito che la media della variabile aleatoria scarto è

nulla: [E(X - E(X)) = 0].

19

2- Sia X una variabile aleatoria discreta; si definisce varianza di X, e si

indica con σ2 (X), la media del quadrato della variabile aleatoria

scarto:

𝜎2 X = E [ ( X − E X )2 ]

Vale a dire:

𝜎2 𝑋 = ( 𝑥𝑖 − 𝐸 𝑋 )2

𝑛

𝑖=1

𝑝𝑖

Le principali proprietà della varianza della variabile X sono le

seguenti:

a) Data X, di media E(X), risulta:

𝜎2 𝑋 = 𝐸 𝑋2 – 𝐸(𝑋)2

b) Se X è una variabile aleatoria, di media E(X), e a e b sono

costanti, risulta:

𝜎2 𝑎𝑋 + 𝑏 = 𝑎2𝜎2 𝑋

2.5. Momenti

Nella pratica non è possibile o non è necessario conoscere l‟intera legge

spaziale della funzione aleatoria, per cui si cerca di individuarne solo alcuni

parametri o momenti. Più precisamente ci si sofferma al calcolo dei primi due

momenti.

- momento del primo ordine

Il momento del primo ordine considerato è il valore atteso (o speranza

matematica). Il valore atteso di una variabile aleatoria Z(x) è definito come:

𝑚 𝑥 = 𝐸 𝑍 𝑥 𝑥𝑜𝑛 𝐸 𝑍 𝑥 < ∞

20

Il valore atteso è una media pesata delle possibile uscite della variabile

aleatoria, con il peso di ogni uscita determinato dalla sua probabilità di

occorrenza. Nel caso di una variabile aleatoria continua Z(u) con possibili

realizzazioni z comprese nell‟intervallo [a,b] e una con f(u; z) la sua p.d.f., la

speranza matematica è espressa come

𝑚 𝑥 = 𝑧.𝑓 𝑥; 𝑧 𝑑𝑧𝑏

𝑎

- momenti del secondo ordine

La varianza di una variabile aleatoria è espressa come il valore atteso degli scarti

quadratici della

media:

𝜎2 𝑥 = 𝑉𝑎𝑟 𝑍 𝑥 = 𝐸 𝑍 𝑥 − 𝐸 𝑍 𝑥 2 = 𝐸{𝑍 𝑥 2 − 𝐸{𝑍 𝑥 }2

Posto che 𝐸{𝑍 𝑥 }2 < ∞

I momenti del secondo ordine considerati sono la covarianza (o

covariogramma), il correlogramma e il variogramma.

Covarianza = 𝐶 𝑥, 𝑥´ = 𝐸 𝑍 𝑥 .𝑍(𝑥2) = 𝐸{𝑍(𝑥) − 𝐸{𝑍(𝑥´)}

Correlogramma = 𝜌 𝑥, 𝑥´ =𝐶(𝑥 ,𝑥 ´)

𝐶 𝑥 ,𝑥 .𝐶(𝑥 ´,𝑥 ´)

Variogramma = 2𝛾 𝑥, 𝑥´ = 𝑉𝑎𝑟{𝑍 𝑥 − 𝑍(𝑥´)}

Osservando la formulazione dei momenti del processo stocastico, ci si può

chiedere come tali parametri possano essere ricavati dai dati disponibili, i quali

rappresentano soltanto alcuni frammenti di un‟unica realizzazione della funzione

aleatoria. Nella pratica la stazionarietà (Cressie, 1993) e richiesta solo per i

momenti di primo e secondo ordine.

21

Tornando al calcolo del momento di primo ordine l‟ipotesi di ergodicità

garantisce che ciò che calcoliamo dai nostri dati sia rappresentativo del processo

stocastico: dall‟analisi di un‟altra realizzazione avremmo ottenuto i medesimi

risultati statistici. L‟ipotesi di stazionarietà invece ci dice che il momento del primo

ordine è stazionario è quindi è lo stesso qualsiasi sia la posizione x considerata.

Ora possiamo utilizzare i dati disposti nello spazio come se fossero diverse

realizzazioni di una singola variabile aleatoria: il calcolo del valore atteso diviene

possibile.

𝑚 = 𝐸{𝑍(𝑥)}

Applicando l‟ipotesi di stazionarietà la formulazione dei momenti del primo

e secondo ordine può essere modificata ad indicare la loro indipendenza dalla

posizione. La funzione covarianza, il variogramma e il correlogramma non sono

più dipendenti dalla posizione assoluta dei punti x e x‟ considerati ma dal vettore

di separazione h.

𝜎2 = 𝑉𝑎𝑟 𝑍 𝑥

Covarianza = 𝐶 ℎ = 𝐸 𝑍 𝑥 .𝑍 𝑥 + ℎ − 𝐸 𝑍 𝑥 .𝐸{𝑍(𝑥 + ℎ)}

Correlogramma = 𝜌 ℎ = 𝐶(ℎ)

𝐶(0)

Variogramma 2𝛾(ℎ) = 𝑉𝑎𝑟 𝑍 𝑥 − 𝑍 𝑥 + ℎ

Si osservi che le ipotesi di stazionarietà ed ergodicità sono relative al

modello e non hai dati ed in quanto tali non sono dimostrabili mediante l‟analisi

dei dati. Nella pratica non viene usata l‟ipotesi di stazionarietà in senso stretto ma

ne vengono utilizzate delle formulazioni meno stringenti. Elenchiamo qui di

seguito le principali definizioni di stazionarietà:

(a) Stazionarietà in senso stretto

E‟ l‟ipotesi più stringente in cui si ammette che l‟intera legge

spaziale sia indipendente dalla posizione considerata.

22

Nella pratica non viene usata l‟ipotesi di stazionarietà in senso

stretto ma ne vengono utilizzate delle formulazioni meno stringenti.

(b)Stazionarietà di secondo ordine

Con questa ipotesi si ammette che i momenti del primo e del

secondo ordine esistano e siano indipendenti dalla posizione considerata.

Più precisamente si dice che una funzione aleatoria Z(x) è stazionaria del

primo ordine se:

𝐸 𝑍 𝒙 = 𝑚 ∀ 𝒙 Є 𝐴

Il valore atteso esiste e non dipende dalla posizione 𝒖 considerata

𝐶 ℎ = 𝐸 𝑍 𝒙 . 𝑍(𝑥 + ℎ) − 𝑚2 ∀ 𝒙 Є 𝐴

Il covariogramma esiste e dipende solo dal vettore 𝒉.

Da queste semplici proprietà ne discendono poi delle altre molto

importanti:

𝐶 0 = 𝑉𝑎𝑟 𝑍 𝒙 ∀ 𝒙 Є 𝐴

𝐶 𝒉 = 𝐶 0 − 𝛾 𝒉

(C) Stazionarietà intrinseca

Per cercare di superare il problema delle misure non stazionarie e

con varianza indefinita, Matheron (1963, 1965) introdusse un‟ipotesi di

lavoro, l‟ipotesi intrinseca, che non prende direttamente in

considerazione la variabile, ma i suoi incrementi: gli incrementi della

funzione stocastica considerata sono stazionari

𝐸 𝑍 𝒙 − 𝑍 𝒙 + 𝒉 = 0

La differenza attesa tra i valori misurati in due siti separati da una

distanza h è funzione solo di h, inoltre essa sarà uguale a zero

𝑉𝑎𝑟 𝑍 𝒙 − 𝑍(𝒙 + 𝒉) = 2𝛾(𝒉)

La varianza delle differenze fra osservazioni separate da una

distanza h è finita e dipende solo dal valore della distanza tra i siti e non

dalla posizione locale.

23

La funzione γ(h) e chiamata semi-variogramma. In questo caso ci si

riferisce ad incrementi di ordine 0 ma l‟ipotesi intrinseca può essere

generalizzata ad incrementi di ordine superiore K, parlando, in questo

caso, di stazionarietà intrinseca di ordine K.

Tale concetto verrà ripreso nel paragrafo riguardante il Kriging con trend;

si ricorda, comunque, che per ordine K si intende il grado della deriva che gli

incrementi considerati riescono a filtrare. Ad esempio gli incrementi di ordine 0

permettono di filtrare una deriva costante, ovvero di grado 0.

2.6. CONTINUITA SPAZIALE

Misure di continuità spaziale quantificano la relazione tra il valore di una

variabile in un punto e il valore della stessa variabile (o di un‟altra variabile

differente) in un altro punto.Il concetto di base della Continuità Spaziale è

semplice: più due valori sono vicini più i valori delle variabili in essi misurati si

somigliano (Issaks e Srivastava, 1989). I più diffusi strumenti di analisi della

continuità spaziale sono i seguenti: h-scatterplot, Variogramma, misure di

covarianza e correlogramma. Un h-scatterplot mostra tutte le possibili coppie di

valori di dati le cui posizioni sono separati da una certa distanza h in una

direzione particolare. La posizione del punto (xi, yi) è indicata come ti, e la

separazione tra due punti ti e tj può essere indicata con hij, (figura 2.3)

Fig. 2.3: Esempio per coppie di valori di dati separati da un certa distanza h

24

La forma della nuvola di punti su un grafico a h-scatterplot ci dice come i

valori dei dati sono continui ad una certa distanza in una direzione particolare. Se

i valori in punti distanti h sono molto simili, le coppie si trovano vicino alla linea x

= y, diventano meno simili, la nuvola di punti diventa più diffusa e dispersa.

Guardando i quattro h-scatterplot (figura 2.4), vediamo che da h = (0,1) ad

h = (0,4), la nube diventa sempre più stesa, e meno simmetrica rispetto alla retta

di 45 gradi.

Anche la direzione considerata è molto importante nella costruzione di un

h-scatterplot. Nella figura (2.4), ci sono alcuni punti che non complotto vicino al

resto della nube. per esempio, in ciascuna delle scatterplots c'è un punto che il

fondo rustico del scatterplot. in ogni caso questo punto comporta la 19 ppm del

valore di V, situato ai margini settentrionali del nostro 10 x 10 m2.

Nella figura (2.4a), questo valore è associato con il valore 77 ppm

immediatamente a sud di contabilità per il punto di insolito che riporta lontano

dalla linea di 45 gradi. Dall'altro h-scatterplot nella stessa figura, il valore di 19

ppm è associato con 91 ppm in valore (2.4b), il valore a 64 ppm in (2.4c), ed i

108 ppm valore in (2.4d).

Issaks e Srivastava (1989) affermano che la statistica di sintesi può

essere influenzata notevolmente da alcuni valori dati aberranti (outlier) ed si deve

spesso verificare come esse possono variare dopo che alcuni outlier dati sono

stati rimossi.

Come abbiamo ricordato prima, una delle caratteristiche essenziali di un h-

scatterplot è la dispersione della nuvola di punti: questa caratteristica viene

riassunta nel coefficiente di correlazione. All‟aumentare della dipersione della

nuvola di punti, il coefficiente di correlazione diminuisce.

Il rapporto tra il coefficiente di correlazione di un scatterplot-h e h è detta

funzione di correlazione o correlogramma. Il coefficiente di correlazione dipende

da h che, essendo un vettore, ha sia una grandezza che una direzione.

25

Fig. 2.4: Esempio di come i dati possono essere accoppiati per ottenere un grafico a h-

scatterplot

Seguendo le normali prassi geostatistiche, useremo:

X: vettore di coordinate spaziali (con componenti x, y o "est" e "nord").

z(x): variabile in questione in funzione della posizione spaziale.

h: vector ritardo che rappresenta la separazione tra due posizionispaziale.

z(x+ h): versione ritardati della variabile in esame.

26

- Funzione di covarianza (covariogramma)

𝐶 ℎ =1

𝑁 ℎ 𝑧(

𝑁 ℎ

∝=1

𝑥∝). 𝑧 𝑥∝ + ℎ −1

𝑁 ℎ 𝑧 𝑥∝ .

1

𝑁 ℎ

𝑁 ℎ

∝=1

𝑧(𝑥∝ + ℎ)

𝑁(ℎ)

∝=1

Con 𝐶 (ℎ) = 𝐶 (−ℎ)

- variogramma

Ricordando l‟espressione della varianza, il variogramma può essere espresso

come:

2𝛾 ℎ = 𝑉𝑎𝑟 𝑍 𝑥 − 𝑍(𝑥 + ℎ) = 𝐸{ 𝑍 𝑥 − 𝑍(𝑥 + ℎ) 2} − 𝐸 𝑍 𝑥 − 𝑍(𝑥 + ℎ) 2

Data l‟ipotesi di stazionarietà, il secondo termine dell‟equazione è nullo e il

variogramma può essere espresso come:

2𝛾 ℎ = 𝐸 𝑍 𝑥 − 𝑍(𝑥 + ℎ) 2

Per cui:

𝛾 ℎ =1

2𝑁(ℎ) 𝑧 𝑥∝ − 𝑧(𝑥∝ − ℎ)

2

𝑁(ℎ)

∝=1

Con 𝛾 ℎ = 𝛾 (−ℎ)

In queste espressioni N(h) indica il numero di coppie di dati sperimentali

separati dal vettore di separazione h. γ(h) e C(h) indicano che ci si sta riferendo

rispettivamente a delle stime del variogramma e covariogramma sperimentali non

coincidenti con la controparte del modello stocastico. Come vedremo, nel

paragrafo seguente, il passaggio tra misure di continuità sperimentali e momenti

del modello stocastico necessita di un ulteriore passo: la modellazione del

variogramma.

Ora, con N (h), che rappresenta il numero di coppie separate da h (più o

meno la tolleranza), possiamo calcolare le statistiche per h, Figura(2.5).

27

Fig. 2.5: (a) la correlazione cofficienint, (b) la covarianza, e (c) il variogramma.

Covarianza:

𝐶 ℎ =1

𝑁(ℎ) 𝑧 𝑥∝ . 𝑧(𝑥∝

𝑁(ℎ)

∝=1

+ ℎ) −𝑚0.𝑚+ℎ

Correlazione:

𝜌 ℎ =𝐶(ℎ)

𝜎0.𝜎+ℎ

Semivarianza:

𝛾 ℎ =1

2𝑁(ℎ) 𝑧 𝑥∝ + ℎ − 𝑧(𝑥∝)

2

𝑁(ℎ)

∝=1

Dove m0 e m+h sono i mezzi dei valori di coda e la testa:

𝑚0 =1

𝑁(ℎ) 𝑧 𝑥∝

𝑁(ℎ)

∝=1

𝑚+ℎ =1

𝑁(ℎ) 𝑧 𝑥∝ + ℎ

𝑁(ℎ)

∝=1

σ0 e σ+h sono le deviazioni standard corrispondente:

28

𝜎0 =1

𝑁(ℎ) 𝑧 𝑥𝛼 − 𝑚0

2

𝑁(ℎ)

∝=1

𝜎+ℎ =1

𝑁(ℎ) 𝑧 𝑥𝛼 + ℎ − 𝑚+ℎ

2

𝑁(ℎ)

∝=1

Il grafico sopra (figura 2.6), mostra tutte e tre le statistiche rispetto al

ritardo medio per le coppie di dati che contribuiscono a ciascun lag. Il più piccolo

ritardo illustrato (lag "zero") comprende sei coppie di dati con un ritardo medio di

351 m.

Fig. 2.6: mostra covarianza, correlazione e semivarianza rispetto all‟effettivo ritardo

medio per le coppie di dati che contribuiscono in ogni punto.

La correlazione contro la distanza h denominata correlogramma e la

semivarianza contro la distanza è il semivariogramma. La covarianza contro la

distanza è il covariogramma. La stima di queste funzioni sulla base di dati

irregolarmente distribuiti può essere molto difficile a causa della necessità di un

sufficiente numero di coppie di dati in lag.

29

Distanze di lag più grandi e/o con tolleranza più grande aumentano le

coppie di dati per lag per la stima; ma riducono la quantità di dettaglio del

semivariogramma (o covarianza o correlogramma).

Il problema è particolarmente difficile per lag piccoli dove solitamente ci

sono poche coppie di dati. Questo va a variare il comportamento della stima

nelle vicinaze all'origine che talvolta è molto importante. Sotto la condizione di

stazionarietà di secondo ordine (media spazialmente costante e varianza), le

funzioni di covarianza, correlogramma e semivariogramma assumono le seguenti

formule:-

𝐶 0 = 𝐶𝑜𝑣 𝑍 𝑥 ,𝑍 𝑥 = 𝑉𝑎𝑟 𝑍 𝑥

𝜌 ℎ =𝐶 ℎ

𝐶 0

𝛾 ℎ = 𝐶 0 − 𝐶(ℎ)

A parole, la covarianza a lag zero è uguale alla varianza globale della

variabile in esame. Il correlogramma è la funzione di covarianza fratto la

varianza, e il semivariogramma è la funzione di covarianza capovolta.

Ciò indica che la covarianza corrispondente è ottenuta dall‟inverso del

variogramm. La figura (2.7) illustra questa idea. A differenza delle analisi di serie

temporali, che preferiscono lavorare con la funzione di covarianza o il

correlogramma, in geostatistica di solito si lavora con il semivariogramma. Ciò è

dovuto al fatto che il semivariogramma, tende a filtrare la varianzione della

variabile rispetto alla media. Inoltre, il semivariogramma può essere applicato in

caso di differenze della variabile, Z (x) - Z (x + h), nel caso di stazionarietà del

secondo ordine. Questa forma di stazionarietà, di cui come l'ipotesi intrinseca, è

un requisito più debole del secondo ordine stazionarietà della stessa variabile,

nel senso che la semivariogramma può essere definito in alcuni casi in cui la

funzione di covarianza non può essere definita.

30

In particolare, il semivariance può continuare ad aumentare con lag

crescente, anziché stabilizzarsi, corrispondente ad una variazione infinita globale.

In questo caso la funzione di covarianza non è definita.

Fig. 2.7: Considerando che il variogramma comincia da zero ed aumenta fino ad un

limite, la covarianza spaziale comincia fuori dalla varianza e diminuisce.(da Raspa,

1995).

2.7.VARIOGRAMMA

È la funzione più comune in geostatistica è:- permette di caratterizzare la

variabilità spaziale diun fenomeno regionalizzato (Raspa, 1995). Il variogramma

o semivariogramma è unafunzione matematica che permette di prevedere come

il valore misurato o previsto di unavariabile regionalizzata cambi in base alla sua

dislocazione nello spazio bi- o tridimensionale.

31

Permette quindi di definire e quantificare la differenza di valore fra punti di

campionamentoseparati da una data distanza h (Hack, 2005). Il

semivariogramma è un modello statisticobivariato (coppie di punti), e la sua

formulazione è connessa con la covarianza e lacorrelazione spaziale. Come

queste ultime infatti, il semivariogramma stima il grado medio varianza tra valori

in funzione della loro distanza nello spazio geografico (Rossi et al., 1992).

Il semivariogramma numericamente consiste nella semivarianza

dell‟incremento [Z(xi+h) - Z(xi)]:

𝛾 ℎ =1

2𝑁(ℎ) 𝑧 𝑥𝑖 − 𝑧(𝑥𝑖 + ℎ)

2

𝑁(ℎ)

𝑖=1

Dove: Z(xi) = valore della variabile Z nel punto xi; Z(xi+h) = valore della

variabile ; Z in un punto distante h rispetto a xi ;n = numero di coppie di punti

separati da una distanza inferiore o uguale ad h. Per descrivere e caratterizzare

un semivariogramma sperimentale, sia esso omnidirezionale o direzionale, è

necessario adattare ad esso un modello matematico. In questo modo viene

creato un semivariogramma modello.

I grafico (2.8) seguente esprime il comportamento spaziale della variabile

casuale, oltre a mostrare dimensione della zona di influenza e gli aspetti

occasionali di anisotropia. Si prevede che le osservazioni (dati, campionamento)

geograficamente vicini gli uni agli altri sono più simili gli uni agli altri da quelli a

distanze maggiori. Pertanto, si prevede che 𝛾 (h) aumenta con la distanza h.

Figura (2.8) mostra un adattamento del variogramma teorico al variogramma

sperimentale con le caratteristiche vicino al ideale.

Il variogramma modello è definito da 3 parametri fondamentali, i cui valori

sono necessari per una corretta applicazione delle procedure di interpolazione

geostatistica ( sill, range e nugget). Si noti che la varianza della differenza

aumenta con la distanza, in modo che il semivariogramma si può pensare come

una funzione di dissomiglianza.

32

Fig. 2.8: Parametri fondamentali della funzione variogramma

L'altezza che raggiunge il semivariogramma si chiama sill. Esso è si è

composto da due parti: una discontinuità all'origine, chiamata nugget, e il sill

parziale. La distanza alla quale il semivariogramma raggiunge il sill è detto range

(Johnson et al., 2001; Rossi et al., 1992; Studiorum, 2006).

I valori del sill, range e, nugget, forniscono importanti indicazioni riguardo

alla struttura dei dati. Sono utilizzati nel kriging per definire i pesi da attribuire alle

singole misure nell‟interpolazione. Nella figura (2.8) sono riportati i parametri che

definiscono la funzione semivariogrammastazionario:

Il nugget (Cο )rappresenta la discontinuità all‟origine, ed è espressione

delle variazionicasuali, della variabilità su piccola scala (inferiore cioè alla

scala del campionamento),e dell‟errore di campionamento.

Il range (a) corrisponde alla distanza dell‟intervallo h in cui il

semivariogramma (in unmodello stazionario) raggiunge un valore costante

detto sill. Intuitivamente a è la massima distanza per cui due punti

possono essere considerati spazialmente correlati.

33

Il sill (C+Cο) rappresenta la varianza massima tra i punti misurati, che per

modelli stazionari corrisponde alla varianza di campionamento (Hack,

2005).

In alcuni casi di studio si osserva che la continuità e la regolarità spaziale

della VR sono responsabili del comportamento del variogramma nell‟origine (cioè

quando h → 0). Si possono individuare tre forme dicomportamento: parabolico,

lineare e discontinuo. Se il variogramma presenta unadiscontinuità all‟origine,

γ(h) non tende a zero quando h tende a zero, pur dovendo sempreessere per

definizione γ(0) = 0.

Dopo diverse iterazioni del variogramma sperimentale in diverse direzioni, è

necessario regolare il modello matematico per rappresentare la variabile più

realisticamente possibile.

Le stime ottenute mediante il kriging saranno quindi più precise e affidabili.

Diversi modelli teorici (riportati in figura 2.9, 2.10 e 2.11) sono utilizzati per

modelizzare il variogramma speriamentale.

34

Fig. 2.9: Modelli principali di variogramma Teorico

35


36


37

2.8. INTERPOLAZIONE SPAZIALE

La mappatura di qualsiasi tipo di variabile richiede un certo numero di

misure (in punti identificati da un sistema di coordinate) e una tecnica di

interpolazione (deterministica o stocastica) per ricostruire i valori assunti dalla

grandezza in oggetto negli intervalli compresi fra i punti noti. Le tecniche

deterministiche stimano, sulla base delle osservazioni, una combinazione lineare

di funzioni che descrive l'andamento medio di una grandezza senza riprodurne i

valori nei punti noti (metodo dei minimi quadrati), o che assume i valori esatti nei

punti noti, fornendo stime poco attendibili nelle regioni comprese fra questi

(interpolazione polinomiale). Fra le tecniche deterministiche possiamo citare la

trendsurface analisys che si occupa della definizione di trends mono-, bi- o

pluridimensionali sulla base di osservazioni regolarmente distribuite nello spazio.

Per una tecnica di interpolazione stocastica quale il kriging, ogni

osservazione viene considerata come una singola realizzazione di una variabile

aleatoria di cui è noto (o ipotizzato) il valore medio in ogni punto (il trend) e le cui

proprietà statistiche sono definite da una funzione detta variogramma. Sulla base

delle osservazioni disponibili vengono stimati tutti i valori desiderati, mentre quelli

noti sono ricostruiti esattamente. Caratteristica del kriging è la possibilità di

disporre, per ogni valore ricostruito, di una stima dell'affidabilità della

ricostruzione. Ciò consente di definire per quest'ultima un intervallo fiduciale

ovvero individuare le aree in cui è necessario aumentare la densità dei

rilevamenti. Inoltre la definizione del modello di interpolazione non avviene in

maniera del tutto univoca (pur esistendo dei criteri guida, il risultato è affidato in

qualche misura all'abilità del modellista ed alla sua esperienza nel campo

specifico di applicazione). L‟interpolazione spaziale è un problema classico di

stima di una funzione F(xp) dove xp è un punto generico del piano, con coordinate

x e y che si vuole stimare a partire da un certo numero di valori noti „m‟ vicini xi:

38

𝑭 𝒙𝒑 = 𝑾𝒊

𝒎

𝒊=𝟏

.𝑭 (𝒙𝒊 )

Il problema che si deve affrontare è la stima dei pesi 𝑾𝒊 di ciascuno dei punti.

Il scopo dell‟interpolazione consiste in: nel trovare un modello plausibile per

esprimere il fatto che punti vicini hanno mediamente valori simili tipi di modelli di

interpolazione; I metodi di interpolazione possono essere suddivisi in due

categorie:

Deterministici: (Poligoni di Thiessen, IDW, Splines) che si basano

sull‟utilizzo di funzioni matematiche.

Stocastici o geostatistici (Kriging): usano la teoria della probabilità con

l‟ applicazione di un modello statistico (semivariogramma) partendo dalla

covarianza tra i punti in funzione della distanza tra gli stessi e scegliendo i

pesi a partire dal grado di somiglianza tra i diversi valori noti della

funzione..

2.9. METODI DETERMINISTICI

I metodi deterministici sono utilizzati per creare mappe della variabile da

analizzare a partire dai valori misurati, basandosi sul grado di similarità o sul

grado di smoothing tra i punti. Tra le tecniche deterministiche si possono

annoverare (Hack, 2005):

• Trend-surface analysis (TSA):- enfatizza la configurazione della

distribuzione spaziale dei dati attraverso lo smoothing delle irregolarità

locali per evidenziare eventuali trend altrimenti celati dal rumore di fondo.

Matematicamente i dati strutturati spazialmente vengono divisi nelle due

componenti locale e globale.

39

• Geometric methods (area interpolation method):- tra queste si possono

citare

Polygonal estimation weighting, cell clustering e Nearest neighbor

estimation.

• Inverse distance weighting (IDW).(point interpolation method):- si basa

sull‟osservazione empirica che il valore della media pesata di un campione

in un determinato punto è proporzionale al peso inverso della distanza dei

campioni che lo circondano.

La tecnica di interpolazione IDW si basa sull‟assunto generale per cui le

misure in punti che sono vicini tendono ad essere più simili di quelle prodotte in

punti più lontani tra loro. Il metodo dell‟IDW ricostruisce la distribuzione spaziale

di una variabile attraverso un modello matematico espresso da una superficie

costituita da punti interpolati. Per stimare il valore della variabile in qualsiasi

punto in cui non sono disponibili misure, il metodo dell‟IDW utilizza i dati

disponibili (misurati) nell‟intorno del punto da interpolare.

Le misure più vicine al punto da stimare eserciteranno una maggiore

influenza sul valore interpolato rispetto a quelli posti a maggiore distanza. Il metodo dell‟IDW assume che ogni misura ha un‟influenza locale che diminuisce

con la distanza. E anche assegna ai dati utilizzati nell‟interpolazione un peso

decrescente allontanandosi dal punto da interpolare.

La formula generale di interpolazione utilizzata nel metodo dell‟IDW è

quella della media pesata:

𝑍 𝑥0 = 𝜆𝑖 𝑍

𝑁

𝑖=1

𝑥𝑖

- 𝑍 𝑥0 è il valore della variabile Z da stimare nel punto (punto previsione);

- i = 1, 2, 3,....., N, è il numero di punti utilizzati nella stima per

interpolazione; si tratta di punti che circondano il punto da stimare e per i

quali sono disponibili misure della variabile Z;

40

- 𝜆𝑖 è il peso assegnato a ciascuna misura spermentale usato

nell‟interpolazione;

- 𝑍 𝑥𝑖 è il valore misurato della variabile Z nel punto 𝑥𝑖 .

Nel metodo dell‟IDW la somma di tutti i pesi applicati alle misure utilizzate per

l‟interpolazione è uguale a 1:

𝜆 = 1

𝑁

𝑖=1

(unto s

2.10. METODI STOCASTICI

Tra le tecniche stocastiche rientra l‟interpolazione Kriging, una tecnica a

media pesata mobile, che per alcuni aspetti può essere ricondotta

all‟interpolazione IDW. La differenza principale risiede nel fatto che oltre a

considerare la distanza dei punti di campionamento, la tecnica Kriging prende in

considerazione anche la disposizione di tali punti nello spazio.

La stima del valore in un punto non campionato viene fatta mediante

l‟utilizzo della funzione variogramma, e la configurazione spaziale dei punti di

campionamento e rispettivi valori misurati attorno al punto da stimare. Tale

tecnica quindi ha come assunto di base che i dati siano parte di un modello

stocastico stazionario, e che i dati siano normalmente distribuiti (Ramilli, 2006).

Talvolta i metodi deterministici (es., IDW) hanno una accuratezza

paragonabile al metodo Kriging, ma sono molto più semplici e pratici e non

presentano gli assunti di base restrittivi delle tecniche stocastiche. Prima di

passare alla modellistica spaziale quindi, l‟analisi dei dati dovrebbe essere

condotta mediante l‟impiego di metodi deterministici (Liebhold et al., 1993; Weisz

et al., 1995).

41

2.11. IL KRIGING

In ambito minerario, ai fini della valutazione della produttività dei

giacimenti, vi era la necessità di una serie di algoritmi di stima che permettessero

di operare con datasets caratterizzati da una geometria di campionamento

irregolare e che tenessero conto del duplice aspetto strutturato - casuale dei

fenomeni analizzati, in questo caso rappresentati dalla distribuzione spaziale di

concentrazioni di metalli pregiati. L‟ingegnere sudafricano Krige pensava a

questo, quando ha costruito il suo algoritmo di stima: di qui il termine Kriging.

Con il termine Kriging si indica una famiglia di stimatori di regressione lineare

non distorti e a minima varianza. La stima 𝒛 ∗ 𝒙 in un punto 𝒙 viene compiuta

utilizzando una combinazione lineare degli n dati noti 𝒛 𝒙∝ , qui concettualizzati

come realizzazioni delle n variabili aleatorie 𝒛 𝒙∝ :

𝑧∗ 𝑥 = 𝜆∝ 𝑥 𝑧(𝑥∝)

𝑛

∝ =1

Dove 𝜆∝ indica il peso attribuito al dato 𝑧(𝑥∝).

La scelta ideale dei pesi sarebbe quella che permette una stima perfetta,

ovvero tale per cui l‟errore di stima- 𝒛 ∗ 𝒙 − 𝒛(𝒙), dove 𝒛 𝒙 è il valore vero ( e

ignoto), sia nulla. Ovviamente ciò non può essere compiuto dato che non

conosciamo 𝒛 𝒙 . Ma, trattando l‟errore di stima in termini di variabili aleatorie,

possiamo cominciare a ragionare sui momenti di primo e secondo ordine

dell‟errore. Dalle stime passiamo allo stimatore di Kriging che, essendo una

combinazione lineare di variabili aleatorie. Esso stesso diviene una variabili

aleatoria:

𝒁∗ 𝒙 = 𝝀∝ 𝒁(𝒙∝)𝒏∝ =𝟏 (k1)

I pesi 𝝀∝ vengono scelti in modo tale che lo stimatore sia ottimale e non

distorto. L‟ottimalità viene raggiunta minimizzando la varianza dell‟errore (Isaaks

e Srivastava, 1989):

42

𝑽𝒂𝒓 𝒁 ∗ 𝒙 − 𝒁 𝒙 = 𝒎𝒊𝒏𝒊𝒎𝒂

la non distorsione si ottiene imponendo la speranza matematica nulla:

𝑬 𝒁 ∗ 𝒙 − 𝒁 𝒙 = 𝟎

Molti descrivono la forma generica dello stimatore del Kriging nella seguente

forma:

𝒁∗ 𝒙 − 𝒎 𝒙 = 𝝀∝ [𝒁 𝒙∝ −

𝒏

∝ =𝟏

𝒎 𝒙∝ ] (k2)

𝒁∗ 𝒙 = 𝝀∝ 𝒁 𝒙∝ +

𝒏

∝ =𝟏

(𝒎 𝒙 − 𝝀∝ 𝒎 𝒙∝ ) (k3)

E ricordandosi del dualismo trend- residuo

𝑹∗ 𝒙 = 𝝀∝ (𝒙)𝑹 𝒙∝

𝒏

∝ =𝟏

(k4)

𝒁∗ 𝒙 = 𝝀∝ 𝒙 𝑹 𝒙∝ + 𝒎(𝒙)

𝒏

∝ =𝟏

(k5)

La (k2) si riduce alla (k1) quando la deriva è ignota, poichè in questi casi

viene filtrata imponendo delle condizioni sui pesi. L‟equazione del kriging

espressa nella forma (k2 o k4) viene utilizzata per la sua generalità e per

evidenziare come la condizione di “non distorsione” sia ottenibile cercando il

filtraggio della deriva (e viceversa) e che questo avvenga ponendo dei vincoli sui

pesi della combinazione lineari. Chiaramente, nella pratica, la forma utilizzata è la

(k1), dato che il più delle volte la media 𝒎(𝒙), essendo ignota, viene filtrata.

Come vedremo in seguito, la varianza dell‟errore, che è una combinazione

lineare di variabili aleatorie, è espressa in termini di covarianza. Di qui

l‟importanza nello studio della struttura di continuità, ovvero del momento del

secondo ordine.

43

2.11.1. KRIGING SEMPLICE

Il kriging semplice è utilizzabile solo in quei casi sia possibile ammettere

stazionarietà del secondo ordine e la deriva costante m sia nota:

𝑚 𝑥 = 𝑚 con 𝑚 nota

Dalla (k3) lo stimatore di Kriging diventa:

𝒁∗ 𝒙 = 𝝀∝ 𝒁 𝒙∝ + 𝟏 − 𝝀∝

𝒏

∝ =𝟏

𝒎

𝒏

∝ =𝟏

In questo caso ho sempre 𝑬 𝒁 ∗ 𝒙 − 𝒁 𝒙 = 𝟎 per le assunzioni sulla

media, infatti:

𝑬 𝒁 ∗ 𝒙 − 𝒁 𝒙 = 𝑬 𝝀∝ 𝒁 𝒙∝ + 𝟏 − 𝝀∝

𝒏

∝ =𝟏

𝒎 − 𝒁 𝒙

𝒏

∝ =𝟏

=

= 𝑬 𝝀∝ 𝑹 𝒙∝ + 𝝀∝

𝒏

∝ =𝟏

𝒎 𝟏 − 𝝀∝

𝒏

∝ =𝟏

𝒎 − 𝑹 𝒙 − 𝒎

𝒏

∝ =𝟏

𝑬 𝑹 ∗ 𝒙 − 𝑹 𝒙

= 𝟎

Quindi l‟errore di stima 𝒁 ∗ 𝒙 − 𝒁 𝒙 può essere espresso come:

𝒁 ∗ 𝒙 − 𝒁 𝒙 = 𝑹 ∗ 𝒙 − 𝑹 𝒙 = 𝝀∝ 𝑹

𝒏

∝ =𝟏

𝒙∝ − 𝑹(𝒙)

Ponendo 𝝀∝ = −𝟏 e 𝒙𝟎 = 𝒙 (cioè corrispondente al punto di stima)

l‟errore di stima può essere scritto come:

𝑹∗ 𝒙 − 𝑹(𝒙) = 𝝀∝ 𝑹 𝒙∝

𝒏

∝ =𝟎

44

Questa notazione è particolarmente conveniente in quanto si può

esprimere la varianza dell‟errore nel modo seguente:

𝑽𝒂𝒓 { 𝝀∝ 𝑹 𝒙∝ } + 𝝀∝

𝒏

𝜷 =𝟏

𝝀𝜷

𝒏

∝ =𝟎

𝑪𝑹 𝒙∝ − 𝒙𝜷

𝒏

∝ =𝟎

Da questa formulazione, si comprende che la 𝐶𝑅 ℎ , covarianza dei

residui, deve essere tale da assicurare che la varianza sia maggiore uguale a

zero, il che implica che deve essere una funzione continua positiva definita. La

minimizzazione della varianza dell‟errore avviene ponendo a zero le n derivate

parziali della varianza rispetto agli n parametri 𝝀∝ :

𝟏

𝟐

𝝏𝑽𝒂𝒓{ 𝝀∝ 𝑹 𝒙∝ }

𝝏𝝀∝ 𝝀𝜷

𝒏

𝜷 =𝟎

𝑪𝑹 𝒙∝ − 𝒙𝜷 = 𝟎

Per ogni ∝ = 𝟏,… . . ,𝒏

Giungendo al ben noto sistema di equazioni normali della regressione

lineari:

𝝀𝜷

𝒏

𝜷 =𝟏

𝑪𝑹 𝒙∝ − 𝒙𝜷 = 𝑪𝑹 𝒙∝ − 𝒙

∝ = 𝟏,… . . ,𝒏

In questo caso, grazie alle assunzioni sulla media, sussiste completa

coincidenza tra covarianza dei residui 𝑪𝑹 𝒉 e quella dei dati grezzi 𝑪 𝒉 . Quindi

il sistema di Kriging può essere scritto come:

𝝀𝜷

𝒏

𝜷 =𝟏

𝑪 𝒙∝ − 𝒙𝜷 = 𝑪 𝒙∝ − 𝒙

∝ = 𝟏,… . . ,𝒏

Quindi il modello di funzione covarianza diviene necessario al fine della

risoluzione sistema di equazioni del Kriging. Una caratteristica del Kriging è nello

stimare il valore della variabile in un punto 𝒙 ne fornisce anche la varianza

dell‟errore di stima.

45

Nel caso del Kriging simplice la varianza minimizzata di stima è:

𝝈𝟐 𝒙 = 𝑪 𝟎 − 𝝀∝

𝒏

∝ =𝟏

𝑪 𝒙∝ − 𝒙

Si osservi che nella gran parte degli algoritmi numerici, i sistemi di Kriging

sono generalmente espressi in termini di covariogramma e non variogramma,

mentre l‟inferenza statistica viene compiuta in termini di variogramma (

Goovaerts, 1997). Poichè nell‟ipotesi di stazionarietà intrinseca il covariogramma

non esiste, in fase di calcolo viene ricavato un pseudo-covariogramma,

utilizzando la seguente formula ( con A costante arbitraria):

𝐶 ℎ = 𝐴 − 𝛾 ℎ ≥ 0,∀ℎ

2.11.2. KRIGING ORDINARIO

In certi casi si può assumere che la media m, entro un sottodominio

spaziale 𝑊 𝑥 del dominio studio A, sia costante ma ignota. In questo caso, non

conoscendola, è necessario filtrarla ponendo una condizione sui pesi:

𝒁 ∗ 𝒙 = 𝝀∝

𝒏

∝ =𝟏

𝒁(𝒙∝) 𝟏 − 𝝀∝

𝒏

∝ =𝟏

𝒎

𝒁 ∗ 𝒙 = 𝝀∝

𝒏

∝ =𝟏

𝒁 𝒙∝ 𝒄𝒐𝒏 𝝀∝

𝒏

∝ =𝟏

= 𝟏

In realtà, la condizione sui pesi poteva essere ricavata esprimendo

l‟algoritma del Kriging partendo dalla (k1) e cercando la non distorsione

dell‟errore:

𝒁 ∗ 𝒙 = 𝝀∝

𝒏

∝ =𝟏

𝒁 𝒙∝

46

= 𝑬{𝒁 ∗ 𝒙 − 𝒁(𝒙)} = 𝑬{ 𝝀∝

𝒏

∝ =𝟏

𝒁 𝒙∝ − 𝒁(𝒙)} =

= 𝑬 𝝀∝

𝒏

∝ =𝟏

𝑹 𝒙∝ − 𝑹 𝒙 + 𝝀∝

𝒏

∝ =𝟏

𝒎−𝒎 =

= 𝝀∝

𝒏

∝ =𝟏

𝒎−𝒎

Quindi solo se la sommatoria dei pesi è uguale ad 1, l‟errore risulta non

distorto. Procedendo in entrambe le maniere, partendo dalla (k3) (per poi arrivare

alla (k1)) o partendo dalla (k1) si ottiene il medesimo risultato: la non distorsione

dell‟errore e il filtraggio della media ignota. Sotto la condizione di non distorsione,

come per il kriging semplice, possiamo descrivere l‟errore di stima

indifferentemente come una combinazione lineare in termini di 𝒁(𝒙) che di 𝑹(𝒙),

per cui la varianza di stima ha la medesima formulazione:

𝑽𝒂𝒓 𝝀∝

𝒏

∝ =𝟎

𝒁 𝒙∝ = 𝑽𝒂𝒓 𝝀∝

𝒏

∝ =𝟎

𝑹 𝒙∝ = 𝝀∝

𝒏

𝜷=𝟎

𝝀𝜷

𝒏

∝=𝟎

𝑪𝑹 𝒙∝ − 𝒙𝜷

Sotto la condizione che:

𝝀∝

𝒏

∝ =𝟏

= 𝟏

O ponendo 𝝀∝ = 𝟏

𝝀∝

𝒏

∝ =𝟎

= 𝟎

In questo caso si può dimostrare che la varianza dell‟errore (Journel, 1989,

p.20), sempre e solo per quelle combinazioni lineari che soddisfano le condizioni

sui pesi, può essere descritta in termini di variogramma:

47

𝑽𝒂𝒓 𝝀∝

𝒏

∝ =𝟎

𝒁 𝒙∝ = 𝑽𝒂𝒓 𝝀∝

𝒏

∝ =𝟎

𝑹 𝒙∝ = 𝝀∝

𝒏

𝜷=𝟎

𝝀𝜷

𝒏

∝=𝟎

𝜸 𝒙∝ − 𝒙𝜷

𝒄𝒐𝒏 𝝀∝

𝒏

∝ =𝟎

= 𝟎

Il variogramma risulta uno strumento utile poichè, a differenza della

covarianza, oltre ad esistere in casi di varianza infinita, ci permette di calcolare la

varianza dell‟errore senza conoscere il valore della media locale m. Infatti, il

variogramma essendo calcolato su incrementi 𝑍 𝑥 − 𝑍 𝑥 + ℎ filtra

automaticamente la media costante entro il sottodominio 𝑊 𝑥 .

2𝛾 = 𝐸 𝑍 𝑥 − 𝑍 𝑥 + ℎ 2 = 𝐸 𝑅 𝑥 − 𝑅 𝑥 + ℎ 2

La minimizzazione della varianza deve avvenire tenendo conto di

condizioni aggiuntive rispetto al Kriging semplice. Affinchè il sistema di equazioni

sia risolvibile si è costretti a ricorrere al metodo dei moltiplicatori di Lagrange µ. I

pesi devono soddisfare il seguente sistema di equazioni affinché la varianza

dell‟errore sia minimizzata:

𝝀𝜷

𝒏

𝜷=𝟏

𝑪𝑹 𝒙∝ − 𝒙𝜷 + 𝑪𝑹 𝒙∝ − 𝒙

∝ = 𝟏,…… ,𝒏

𝝀𝜷

𝒏

∝=𝟏

= 𝟏

Di nuovo, poichè 𝑚 𝑥 è assunta costante all‟interno del dominio 𝑊 𝑥

abbiamo coincidenza tra 𝐶𝑅 ℎ e 𝐶 ℎ . Inoltre il sistema di Kriging può essere

espresso in termini di variogramma (Goovaerts, 1997):

48

𝝀𝜷

𝒏

𝜷=𝟏

𝜸 𝒙∝ − 𝒙𝜷 − µ 𝒙 = 𝜸 𝒙∝ − 𝒙

∝ = 𝟏,…… ,𝒏

𝝀∝

𝒏

∝=𝟏

= 𝟏

Nel caso del Kriging ordinario la varianza minimizzata di stima è la seguente:

𝝈𝟐 𝒙 = 𝑪 𝟎 − 𝝀∝

𝒏

∝=𝟏

𝑪 𝒙∝ − 𝒙 − µ 𝒙

2.12. COKRIGING

In alcuni casi si osserva l'esistenza di una certa correlazione spaziale tra

diverse proprietà distribuite nello spazio. Ad esempio l'impedenza acustica

ricavata da un rilievo sismico è correlata alla porosità della roccia, un certo

elemento chimico può essere correlato con altri elementi chimici.

Ancora, il tenore di un certo metallo nel suolo può essere legato al tipo di

subtrato roccioso o al tipo di processi antropici che agiscono in una certa area.

Risulta quindi ragionevole eseguire le stime della proprietà di interesse, non solo

utilizzando i punti in cui tale proprietà è misurata ma anche l'informazione

derivante da misurazioni di altre proprietà, ad essa correlate.

49

L'informazione derivante da misurazioni dirette della proprietà di interesse

viene chiamata informazione primaria. L'informazione derivante da misurazioni di

proprietà correlate alla proprietà di interesse, viene denominata informazione

secondaria.

Ogni proprietà secondaria potrà essere descritta in termini di realizzazioni

di una variabile aleatoria, che chiameremo variabile aleatoria secondaria. Avendo

N variabili aleatorie secondarie, l'algoritmo di stima, che prende il nome di

Cokriging ordinario, può essere espresso come:

𝒁𝟏∗ 𝒙 = 𝝀∝𝒊

𝒏

∝𝒊 =𝟏

(𝒙)

𝑵𝒗

𝒊=𝟏

𝒁𝒊 𝒙∝𝒊

Con

𝝀𝜷𝟏

𝒏𝟏

𝜷𝟏=𝟏

= 𝟏

𝝀𝜷𝒊

𝒏𝒊

𝜷𝒊=𝟏

= 𝟎 𝒊 = 𝟐,…… ,𝑵𝒗

Come prima, la varianza dell‟errore viene espressa come una

combinazione lineare di variabile aleatorie e risulta determinata da una

combinazione lineare di covarianze dei residui:

𝑽𝒂𝒓 𝝀∝𝒊

𝒏𝒊

∝𝒊 =𝟏

𝑵𝒗

𝒊 =𝟏

𝒁𝒊 𝒙∝𝒊 = 𝝀∝

𝒏

𝜷=𝟎

𝝀𝜷

𝒏

∝=𝟎

𝑪𝒊𝑹 𝒙∝ − 𝒙𝜷

𝑵𝒗

𝒊 =𝟏

Tale varianza viene minimizzata sotto le condizioni di non distorsione

ottenendo il sistema di cokriging ordinario:

50

𝝀𝜷𝒋

𝑵𝒗

𝜷𝒋 =𝟏

𝑪𝒊𝒋

𝑵𝒗

𝑱=𝟏

𝒙∝𝒊 − 𝒙𝜷𝒋 + µ𝒊 𝒙 = 𝑪𝒊𝟏 𝒙∝ − 𝒙

∝ = 𝟏,… ,𝒏∝𝒊 𝒊 = 𝟏,… . ,𝑵𝒗

𝝀𝜷𝟏

𝒏𝟏

𝜷𝟏 =𝟏

= 𝟏

𝝀𝜷𝒊 = 𝟎

𝒏𝒊

𝜷𝒊 =𝟏

𝒊 = 𝟐,… . ,𝑵𝒗

Dalle espressioni precedenti si osserva come l'algoritmo sia identico a

quello del kriging ordinario. Lo stimatore è sempre dato da una combinazione

lineare di variabili aleatorie tale per cui l'errore sia minimo e non distorto; con la

differenza che ora le variabili aleatorie rappresentano più proprietà distribuite

nello spazio. Dall'espressione della varianza dell'errore si può osservare come,

ora, sia necessario conoscere non solo i covariogrammi per ognuna delle variabili

aleatorie ma anche i cross-covariogrammi corrispondenti a tutte le combinazioni

di proprietà.

Chiaramente, la combinazione di modelli di continuità diretti ed incrociati

deve essere tale da garantire la non negatività della varianza (Goovaert, 1997),

rendendo difficile la procedura di inferenza dei modelli di continuità diretta ed

incrociata (modello di coregionalizzazione).

Le difficoltà operative rendono il cokriging uno strumento realmente

vantaggioso qualora la proprietà primaria risulti sotto campionata rispetto alle

proprietà secondarie e quando le strutture di continuità delle diverse proprietà

manifestino caratteristiche, nel senso del tipo di modello e dei parametri che lo

definiscono, notevolmente differenti.

51

Non ci si vuole soffermare oltre su questo algoritmo, rimandando a

Goovaerts, (1997) e Isaaks e Srivastava, (1989) per una più attenta discussione

delle sue caratteristiche e varianti. Avremo comunque modo di soffermarci su

una semplificazione del cokriging che è alla base del trattamento

dell'informazione soft.

2.13. KRIGING INDICATORE

Originariamente il kriging indicatore (IK) è una tecnica che nasce per

operare con variabili continue (Journal, 1983a), ma si può facilmente modificare

questa tecnica per farla operare con variabili categoriche ( Bierkens & Burrough,

1993b). Quando usiamo il kriging per studiare fenomeni ad alta variabilità si

presentano problemi. Prendiamo per esempio un campo aleatorio {𝒁 𝒙 , 𝒙 ∈

ℝ𝟐} intrisecamente stazionario su ℝ𝟐e 𝑵 osservazioni del campo {𝒁 𝒙𝒊 =

𝒙𝒊, 𝒊 = 𝟏, . . ., 𝑵}. Se si suppone che sia presente una alta variabilità tra le

osservazioni, allora il semivariogramma empirico non si presta ad alcuna

possibile interpretazione.

Per risolvere il problema di dati molto variabili possiamo pensare di seguire

due strade:

- La prima consiste nell‟eliminare gli outlier presenti nei dati, ciò può essere

fatto tenendo conto di un modello aleatorio della distribuzione dei dati,

spesso legata all‟ipotesi di normalità, questo metodo non è quasi mai la

miglior strada. Infatti spesso accade che siano proprio gli outlier a

contenere le informazioni più importanti e scartarli porterebbe a produrre

stime del tutto inesatte.

- La seconda strada a cui si potrebbe pensare è invece quella di smussare i

dati tramite funzioni di smoothing non lineari ( p.es. 𝑥).

52

Questo tipo di trasformazioni però richiede l‟utilizzo di tecniche di stima

non lineari come il Disjunctive kriging (DK) e il kriging Gaussiano

multivariato (MG) (Journel & Huijbregts, 1978). ln 𝑥 ,

Chiamiamo 𝐃 ⊂ ℝ𝟐 la regione che contiene tutte la 𝑵 osservazioni

disponibili, allora per ogni 𝐱 ⊂ 𝑫 e perso un valore della caratteristica cercata

del campo come soglia 𝐳, definiamo la funzione indicatore come:

𝒊 𝒙; 𝒛 =

𝟏 𝒔𝒆 𝒛 𝒙 ≤ 𝒛

0 𝑎𝑙𝑡𝑟𝑖𝑚𝑒𝑛𝑡𝑖

Lo scopo di questa tecnica è quello di fornire una predizione di una

caratteristica del campo 𝒁(𝒙) in un sito 𝒙𝟎 non visitato, date le 𝑵 osservazioni

fatte su 𝑫 a nostra disposizione. Spesso però sarà di interesse fare questa

predizione non su 𝑫, ma su un suo sottoinsieme 𝑨, che contiene un numero

𝑴 < 𝑁 di osservazioni.

Il Kriging indicatore è un‟ulteriore tipologia di Kriging. Esso non utilizza i

valori propri della variabile da stimare, ma sfrutta delle nuove variabili indicatrici

per ottenere una predizione sui valori. A differenza dei kriging semplice e

ordinario, il kriging indicatore è un tipo di kriging non parametrico e, utilizzando la

funzione di distribuzione cumulata( C.D.F.)

ricostruita tramite delle soglie, per stimare un valore in un punto ignoto.

Il primo passo è quello di stabilire le soglie che verranno poi utilizzate per

creare le C.D.F., in base alla distribuzione dei dati ed a una prima analisi

statistica. Per ogni soglia, la variabile viene trasformata in un valore 0 o 1: 0 nel

caso in cui il valore superi la soglia, 1 nel caso di non superamento.

Il kriging indicatore non stima la variabile indicatrice 𝑰(𝒙𝟎; 𝒛𝒌), bensi esso

fornisce un metodo di stima ai minimi quadrati della C.D.F. alla soglia 𝒛𝒌:

53

Î(𝒙𝟎; 𝒛𝒌) = E{𝒙𝟎; 𝒛𝒌| 𝒏 ]}

Con 𝒏 che rappresenta il condizionamento causato dalla vicinanza di

altri valori nel punto 𝒙𝟎.

Per come viene costruita la funzione aleatoria. indicatrice 𝑰(𝒙; 𝒛), il valore

atteso di 𝑰(𝒙; 𝒛) corrisponde alla C.D.F. della funzione aleatoria 𝒛(𝒙) ; infatti:

𝐸 𝐼 𝑥; 𝑧 = 1 ∗ 𝑃𝑟𝑜𝑏 𝑍 𝑥 ≤ 𝑧 + 0 ∗ 𝑃𝑟𝑜𝑏 𝑍 𝑥 ≤ 𝑧

𝐸 𝐼 𝑥; 𝑧 = 𝑃𝑟𝑜𝑏 𝑍 𝑥 ≤ 𝑧 = 𝐹 𝑧

Perciò la stima della variabile indicatrice 𝑰(𝒙; 𝒛), alla soglia 𝒛, diventerà:

Î 𝒙𝟎; 𝒛 = 𝑬 𝒙𝟎; 𝒛 𝒏 } = 𝑷𝒓𝒐𝒃 { 𝒁 𝒙𝟎 ≤ 𝒛|(𝒏)}

Î 𝒙𝟎; 𝒛 = 𝒘 𝒙𝒊; 𝒛 𝑰 𝒙𝒊; 𝒛

𝒏

𝒊=𝟏

− 𝟏 − 𝒘 𝒙𝒊; 𝒛

𝒏

𝒊=𝟏

𝑭 𝒛

In questo caso gli 𝒘 𝒙𝒊; 𝒛 sono i pesi attribuiti alla soglia 𝒛.

In seguito si passa all‟analisi variografica dei dati delle differenti soglie e si

crea per tutte un modello di variogramma. I vari modelli di variogrammi vengono

utilizzati per ricostruire le C.D.F. in ogni punto da stimare. Avendo ricostruito le

C.D.F. sarà possibile avere, per ogni posizione dell‟area studiata, il valore atteso

o la probabilità di superamento di una data soglia.

54

2.17. CROSS VALIDAZIONE

Dopo aver effettuato delle stime servono dei metodi per verificarne l‟affidabilità,

per confrontare i risultati derivanti dall‟applicazione di metodi differenti e per

valutare l‟effetto nell‟uso di valori diversi dei parametri principali nei vari metodi.

Un metodo molto utilizzato è la cross validation, o validazione incrociata, che usa

i valori misurati come informazione di riferimento rispetto ai valori stimati. Per

implementare un confronto robusto, si rimuovono, uno alla volta, i punti misurati

effettuando la stima in quel punto con tutti gli altri rimanenti. Iterando la

procedura per tutti i valori noti si può confrontare il risultato della stima con le

informazioni misurate. Un metodo molto semplice è lo scatter plot dei punti noti e

delle loro stime (figura 2.12).

Fig 2.12. Cross validation: scatter plot dei punti noti e delle loro stime (a sinistra) e

istogramma dell‟errore standardizzato (a destra)

55

Più i punti sono vicini alla bisettrice e più la stima è accurata. In generale

la dispersione dei valori stimati è minore di quella dei valori reali, perché la loro

varianza è minimizzata dalla procedura di kriging.

L‟istogramma dell‟errore standardizzato:

𝒁 𝒙 − 𝒁∗(𝒙)

𝝈

è uno strumento molto importante per quantificare graficamente la

dispersione (varianza) delle stime (figura 2.12).

Un altro parametro molto utile per quantificare l‟affidabilità della procedura

di stima è la media del quadrato dell‟errore standardizzato:

𝑬{[𝒁 𝒙 − 𝒁∗ 𝒙

𝝈]𝟐}

che dovrebbe essere uguale a uno per una stima perfetta. Questo

parametro è spesso utilizzato per confrontare risultati ottenuti

56

CAPITOLO 3

IL SOFTWARE R

3.1. INTRODUZIONE

Esiste in circolazione una gamma assai vasta di software specializzati

nell'analisi statistica dei dati, basta ricordare SAS, SPSS, STATA,

STATGRAPHICS PLUS, SHAZAM, S PLUS, MINITAB, GAUSS, etc., solo per

citarne alcuni che sono in commercio. Sono prodotti che costituiscono senz'altro

un fondamentale ed insostituibile ausilio per il lavoro dello statistico, tuttavia molti

di questi programmi sono anche alquanto costosi ed è consentito l'uso su licenza

da parte del produttore.

Da alcuni anni a questa parte, soprattutto in ambito universitario (ma non

solo) si sta sempre più diffondendo un nuovo software che merita di sicuro una

debita trattazione e l'interessamento da parte degli statistici e di coloro che si

occupano di analisi dei dati e che costituisce anche una valida alternativa ai

software sopraccitati. Ci stiamo riferendo al software R. R è un ambiente

statistico. Esso è il frutto del lavoro collettivo svolto da un gruppo, sempre più

numeroso, di ricercatori in campo statistico ed informatico a livello mondiale. R

più che un software statistico può essere definito come ambiente, costituito da

una varietà di strumenti, orientato alla gestione, all'analisi dei dati e alla

produzione di grafici, basato sul linguaggio R creato da AT&T Bell Laboratories,

ambiente dal quale è nato un altro software commerciale più noto, S-Plus,

prodotto prima da MathSoft e ora da Insightful.

57

R, a differenza di quest'ultimo è un GNU-Software, ossia è disponibile

gratuitamente sotto i vincoli della GPL (General Public License)3: chiunque può

accedere al codice sorgente e modificarlo, migliorarlo, integrarlo e mettere a

disposizione della comunità scientifica il proprio lavoro nella logica e nella

filosofia del software opensource. La versione iniziale di R venne sviluppata nel

1996 per il sistema operativo MacOS da Ross Ihaka e Robert Gentelman del

dipartimento di Statistica dell'Università di Auckland in Nuova Zelanda4.

Successivamente altri ricercatori iniziarono ad aggiungersi e a fornire il loro

contributo alla scrittura del codice sorgente e al miglioramento dell'applicazione,

sviluppata in seguito anche per l'ambiente UNIX. Nel 1997 è nato l'R

Development Core Team, gruppo formato da statistici di tutto il mondo che si

occupano dello sviluppo e della distribuzione del programma e, ancora più

recentemente, nel 2003 è stata costituita dagli stessi membri dell'R

Develpoment Core Team la R Foundation for Statistical Computing.

Si tratta di una organizzazione non profit che lavorano per il pubblico

interesse allo scopo di promuovere lo sviluppo (attuale e futuro) e la diffusione

del software, di fornire un punto di riferimento per persone, istituzioni ed imprese

che desiderano sostenere il progetto e interagire con la community che ruota

intorno ad R, di gestire e tutelare il copyright di R e della relativa

documentazione.

3.2. CHE COS’è R?

Le sue caratteristiche principali possono essere così riassunte:

a) semplicità nella gestione e manipolazione dei dati;

b) disponibilità di una serie di strumenti per calcoli su vettori, matrici ed

altre operazioni complesse;

58

c) accesso ad un vasto insieme di strumenti integrati per l'analisi statistica;

d) produzione di numerose potenzialità grafiche particolarmente flessibili;

e) possibilità di adoperare un vero e proprio linguaggio di programmazione

orientato ad oggetti che consente l'uso di strutture condizionali e cicliche,

nonché di funzioni create dall'utente.

R si presenta come uno strumento completo e ben si adatta a molteplici

esigenze: è il software maggiormente utilizzato a livello mondiale dai ricercatori in

campo statistico. Vengono ora analizzati i suoi punti di forza ed anche alcune sue

debolezze e/o mancanze.

3.3. LE FUNZIONALITÀ DELL’AMBIENTE R

L'ambiente R è basato sul concetto di "package" tradotto di solito in

italiano con pacchetto. Un package è un insieme di strumenti che svolgono

determinate funzioni, ma può anche contenere solo dati.

Attualmente, per R, è disponibile una vasta gamma di packages

(scaricabili, per il sistema operativo Windows, dall'URL

http://cran.at.rproject.org/bin/windows/contrib) utilizzabili per la risoluzione di

specifici problemi o per analisi statistiche molto particolareggiate. Il cuore di R è

rappresentato dal modulo base (che offre gli strumenti fondamentali per l'analisi

statistica) e attorno a questo modulo "ruotano" una serie di altre librerie

addizionali, alcune delle quali sono già comprese nel programma R al momento

in cui lo si installa, mentre altre librerie ancora, in relazione alle esigenze e

necessità, possono essere aggiunte e installate dall'utente dopo averle scaricate

dal sito succitato. R è un ambiente computazionale molto dinamico ed in continua

evoluzione e frequentemente nuovi ulteriori packages sono realizzati e messi a

disposizione dei ricercatori.

http://cran.at.r/

59

3.4. ALCUNI PACKAGES DI R

In questo paragrafo si vogliono presentare brevemente alcuni packages di

R di particolare interesse per la statistica ambientale:-

GeoR: fornisce una serie di funzioni per l‟analisi di dati geostatistici

GeoRglm: modelli spaziali lineari generalizzati

Gstat: modelli geostatistici univariati e multivariati, previsioni e simulazioni

GRASS:interfaccia tra GRASS 5.0 geographical information system ed il

software R

Spatial: funzioni per kriging e point pattern analysis

Spatstat: spatial Point Pattern analysis

Splancs: Spatial and Space-Time Point Pattern Analysis

Sgeostat: un “object-oriented framework” per i modelli geostatistici

Vardiag: un package per la diagnostica dei variogrammi

3.5. PACKAGE GSTAT

La geostatistica non è un nuovo soggetto nella comunità R, e sono quindi

disponiibili diversi pacchetti con diversi funzionalità ed obiettivi. Alcuni di questi

sono stati sviluppati per scopi didattici, mentre altri hanno caratteristiche molto

più avanzate, tuttavia nella maggior parte di cosi mancano un certo numero di

caratteristiche che sono comunemente utilizzate in geostatistica (Isaaks e

Srivastava, 1989), in particolare il kriging, la modellazione dei variogrammi

multivariabili e il kriging con simulazioni.

60

Alcuni ricercatori si proposero di creare un nuovo software per ovviare agli

inconvenienti dei package precedenti. Il software gstat si presenta come un

programma per la modellizzazione, predizione e simulazione di dati geostatistici

in una, due o tre dimensioni. Gstat `e un programma per la Geostatistica

multivariata, che comprende previsione e simulazione delle variabili.

La prima versione del software fu terminata nel 1992 e già essa offriva

un‟interfaccia interattiva per modellare variogrammi. Nella sua forma originale,

gstat necessitava dell‟ eseguibile e del programma Gnuplot per la visualizzazione

dei grafici. Dal 2003 gstat può essere utilizzato come pacchetto per alcuni soft-

ware statistici molto conosciuti come S-Plus e R.

Le principali funzioni di gstat sono: la stima del variogramma relativo e

incrociato; l‟adattamento di un modello al variogramma stimato interattivamente;

l‟applicazione delle metodologie di stima del kriging ordinario, semplice e

universale in una variabile e per variabili multiple; l‟applicazione del metodo del

cokriging e l‟esecuzione dei processi di simulazione cond

Università degli studi di Padova - Metodologie geostatistiche...

Documents

Transcript of Università degli studi di Padova - Metodologie geostatistiche...