Geostatistica Di Base

DISPENSE DI GEOSTATISTICA APPLICATA

G. Raspa

Capitolo 3 – Geostatistica di base1

1

Il presente capitolo costituisce le dispense del modulo di Geostatistica del corso integrato “Calcolo delle Probabilità e Geostatistica” impartito nel corso di laurea in Ingegneria per l’Ambiente e il Territorio (Laurea di 1° livello) della Facoltà di Ingegneria dell’Università di Roma “La Sapienza”.

3.

GEOSTATISTICA DI BASE .................................................................................3 3.1 Introduzione alla Geostatistica ............................................................................3 3.2 L’approccio probabilistico ..................................................................................9 3.3 Modelli di Base.................................................................................................11 3.3.1 Modelli Stazionari......................................................................................12 3.3.2 Modelli non Stazionari...............................................................................16 3.3.3 Considerazioni sulla Scelta dei Modelli......................................................24 3.4 Analisi ed interpretazione dei variogrammi sperimentali ..................................25 3.4.1 Calcolo dei variogrammi sperimentali .........................................................25 3.4.2 Identificazione del modello di FA ..............................................................28 3.4.3 Comportamento nell’origine del variogramma ...........................................37 3.4.4 Andamento e modelli del variogramma elementare ....................................41 3.4.5 Variogrammi a strutture annidate ...............................................................49 3.4.6 Anisotropie nei variogrammi......................................................................54 3.5 Aggiustamento di un variogramma modello......................................................57 3.6 Regolarizzazione di una FA ..............................................................................61 3.6.1 Calcolo delle funzioni γ ............................................................................66 3.6.2 Regolarizzazione di un effetto pepita..........................................................67 3.7 Varianza di dispersione.....................................................................................69 3.8 Il kriging ordinario............................................................................................75 3.8.1 Gli stimatori lineari .....................................................................................76 3.8.2 Correttezza della stima................................................................................79 3.8.3 Precisione della stima..................................................................................79 3.8.4 Stimatori tradizionali...................................................................................80 3.8.5 Il Kriging Ordinario ....................................................................................82 3.9 Il kriging semplice ............................................................................................94

2

3.

GEOSTATISTICA DI BASE

3.1 Introduzione alla Geostatistica La Geostatistica studia i fenomeni naturali che si sviluppano su base spaziale a partire dalle 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. I metodi della Geostatistica sono applicabili in tutti quei settori delle scienze applicate 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 minerarie, l’idrologia, l’idrogeologia, la scienza dei suoli, l’agronomia, la geotecnica, la geofisica, il telerilevamento, la climatologia, la meteorologia, l' oceanografia, le scienze forestali, la zoologia, l’epidemiologia, l’igiene ambientale. Si consideri un fenomeno spaziale, per es. l’inquinamento di un sito da metalli pesanti. Indicando con z la concentrazione di uno di essi e con x il generico punto di coordinate (xu, xv)1 del campo di indagine, z(x) è una funzione numerica che rappresenta la concentrazione dell' inquinante nei punti del sito. In questo quadro si possono dare le seguenti definizioni: Variabile regionalizzata (VR) Si intende la funzione matematica z(x) sopra introdotta. Il termine regionalizzata specifica che si tratta di una funzione numerica il cui valore dipende dalla localizzazione, espressa normalmente dalle coordinate spaziali, e che si presenta strutturata spazialmente. Campo: E'il dominio nel quale la variabile z è suscettibile di assumere determinati valori e, all' interno del quale, se ne studia la variabilità. Coincide con lo spazio di osservazione (o di indagine) del fenomeno in esame. Supporto: E'l' entità geometrica sulla quale la variabile z è definita od anche misurata; essa è caratterizzata dalle sue dimensioni e dalla sua forma. Quando le dimensioni sono molto piccole (rispetto alla scala del lavoro) il supporto può considerarsi puntuale (per es. un 1 Nel seguito del testo la posizione di un punto nello spazio di lavoro sarà sinteticamente indicata con x, sottointendendo che si tratta di un punto di coordinate xu , xv, essendo u e v i due assi coordinati del riferimento.

3

campione areale di suolo di qualche decina di dm2 può essere considerato puntuale rispetto ad una distanza tra campioni successivi di alcune decine di m). Non è così invece per un pixel di una scena Landsat TM, che ha dimensioni 30 x 30 m, in un’area di indagine di alcuni Km2. Il concetto di supporto e le sue implicazioni giuocano un ruolo importante nella teoria e nelle applicazioni geostatistiche. Data una variabile regionalizzata riferita ad un determinato supporto, si ha che, cambiando la forma e le dimensioni di esso, si ottiene una variabile regionalizzata diversa dalla prima, ma non senza analogia con essa. Si osservi, a titolo di esempio, la figura 3.1. In essa è riportato l’andamento della conducibilità idraulica misurata su campioni successivi di terreno di dimensioni 10 cm; in tratto più grosso è riportato l’andamento dei valori mediati su un metro (su 10 campioni). Nel primo caso la VR ha un supporto 10 cm e nel secondo ha come supporto 1m. Si tratta di due variabili che presentano un andamento molto diverso tra di loro.

Fig. 3.1 - Andamento della conducibilità idraulica riferita a due diversi supporti (da “Stochastic models of fluid flow in heterogeneous media” di Leslie Smith,Proceedings of the workshop on soil spatial variability of the ISSS and the SSSA, Las Vegas, 1984).

In base alla definizione data, una VR è una variabile puramente deterministica; da un punto di vista matematico è una funzione z(x) che assume in ogni punto dello spazio un determinato valore numerico. Si osservi per es. la figura 3.2 che rappresenta l’andamento, lungo una determinata direzione, della tensione acqua-suolo (tas) di un terreno ad un giorno dall’irrigazione. Si può osservare, da una parte, un andamento irregolare alla piccola scala che non incoraggia uno studio matematico diretto, e dall’altra, una variabilità strutturata, cioè una variabilità che sembra ubbidire a delle regole. Vi si notano infatti tratti con elevati valori della tas e tratti con medi e bassi valori. Pertanto, campioni prelevati in vicinanza dei tratti ad alto valore avranno una elevata probabilità di avere un tas elevata, mentre vi avranno una media e bassa probabilità i campioni prelevati negli altri tratti.

4

Fig. 3.2 - Andamento di una variabile spaziale rappresentante la tensione acqua-suolo, lungo una direzione, di un terreno ad un giorno dall’irrigazione (da “Spatial variability of soil-water properties in irrigated soils”, di P.J.Wierenga, Proceedings of the workshop on soil spatial variability of the ISSS and the SSSA - Las Vegas, USA 1984) .

Un approccio corretto allo studio dei fenomeni spaziali deve considerare entrambi gli aspetti della variabilità e fornire degli strumenti operativi alla risoluzione dei problemi. Un tale approccio è quello probabilistico, cioè basato sulle Funzioni Aleatorie. Prima di passare ad illustrare l’approccio probabilistico vediamo come è possibile, facendo uso dei concetti di varianza, covarianza e coefficiente di correlazione empirici, caratterizzare intuitivamente la variabilità spaziale di una VR. Sia z(x) una VR, di supporto puntuale, definita in un’area S avente un’estensione di alcune centinaia di metri (fig. 3.3) e nota in tutti i punti. Si consideri una coppia di punti distanti h(vettore) di posizione x e x+ h, con |h| = h1 piccolo rispetto alle dimensioni di S, per es. un metro (fig. 3.3 A1). In corrispondenza dei due punti la VR assume i valori z(x) e z(x+h). Se, in maniera analoga, si fanno assumere al primo punto della coppia n posizioni xi ( i =1,n ) all’interno di S in maniera che anche il secondo punto stia all’interno del campo, si ottengono altrettante coppie di valori: z(xi), z(xi+h) che graficamente possono essere rappresentati da una nuvola di correlazione in cui gli assi esprimono i valori del primo e del secondo punto (fig. 3.3 A2).

5

Fig. 3.3 - Descrizione empirica della variabilità spaziale: covarianza delle coppie

Essendo |h| piccolo saranno verosimilmente anche piccole le differenze tra le coppie di valori, per cui la nuvola sarà poco dispersa attorno alla retta a 45°, tanto più o tanto meno a seconda della differenza tra i valori di ciascuna coppia. Per una distanza tra le coppie pari a zero, le differenze sono tutte nulle ed i punti della nuvola si collocano esattamente sulla retta a 45°. Quantitativamente l’entità del legame tra le coppie di valori z(xi) e z(xi+h) è misurata dalla covarianza empirica:

Sx , x + h =

n i=1

Z ( xi ) Z ( xi + h) / n −

n i =1

Z ( xi)

n i =1

6

Z ( xi + h) / n / n

od, ancora meglio, dal coefficiente di correlazione empirico:

ρx , x + h =

Sx, x + h Sx ⋅ Sx + h

dove Sx e Sx+h sono le deviazioni standard empiriche dei valori rispettivamente del primo e del secondo punto:

S2x =

n

Z2 ( xi) / n −

i =1

n

Z ( xi )

i =1

n

Z ( xi ) / n / n

i =1

e S2x + h =

n i =1

Z 2 ( xi + h) / n −

n

Z ( x i + h)

i =1

n

Z ( x i + h) / n / n

i =1

Se ora si considera un valore di |h| = h2 più grande di un metro, per es. 5 m (fig. 3.3 B1) è verosimile pensare di ottenere una nuvola di punti (fig. 3.3 B2) più dispersa di quella ottenuta per la distanza di 1 m. La covarianza ed il coefficiente di correlazione saranno più piccoli dei corrispondenti valori . Se ora consideriamo valori di |h| sempre crescenti si ha che la covarianza ed il coefficiente di correlazione decresceranno, potendosi anche annullare se le coppie di valori risultano indipendenti. In tal caso la nuvola di correlazione ha una dispersione uniforme nel piano z(xi) e z(xi+h) (fig. 3.3 C2). Appare pertanto immediato graficare il coefficiente di correlazione (o la covarianza) in funzione della distanza |h| e considerare questa curva come una forma quantitativa che esprime, limitatamente alla direzione in esame, la variabilità spaziale di z(x) (fig. 3.4). Per |h| = 0 le due variabili z(xi) e z(xi+h) coincidono e di conseguenza il coefficiente di correlazione ρ è uguale a uno. Esaminiamo ancora un’altra forma intuitiva di caratterizzare la variabilità spaziale. Sempre con riferimento alle situazioni descritte nella fig. 3.3, consideriamo per ogni coppia di punti xi e xi+h le differenze:

ε( xi , xi + h) = Z ( xi + h) − Z ( x i)

( i = 1, n )

7

r (h) 1

Fig. 3.4 - Correlazione della coppia in funzione della distanza.

0.5

0

h1

h2

h3

h

Queste, per ogni situazione, presentano istogrammi caratterizzati da dispersioni che aumentano all’aumentare di |h|. La fig. 3.5 mostra in termini di varianza degli incrementi ciò che nelle fig. 3.4 era mostrato in termini di coefficiente di correlazione.

Fig. 3.5 - Descrizione empirica della variabilità spaziale: dispersione degli incrementi.

8

Come è noto, la dispersione dei valori attorno alla media è misurata dalla varianza. Pertanto la curva che esprime la varianza degli incrementi in funzione della distanza è un’altra forma per caratterizzare la variabilità spaziale di un parametro regionalizzato. Quando |h| = 0 la varianza è zero e generalmente aumenta all’aumentare di |h| (fig. 3.6).

g (h) S23

Fig. 3.6 - Dispersione degli incrementi in funzione della distanza.

S 22

S 21

h1

h2

h3

h

Il trasferimento di questi concetti di carattere empirico nel quadro di un formalismo teorico esigente costituisce l’approccio probabilistico della Geostatistica allo studio dei fenomeni naturali a carattere spaziale.

3.2 L’approccio probabilistico

Prendiamo in esame un fenomeno spaziale, per es. l’andamento ad un certo istante ed all’interno di un’area S, della tavola d’acqua di un acquifero sedimentario ed indichiamo con z la sua profondità (rispetto al piano di campagna). La funzione matematica, che esprime l’andamento di z in funzione della posizione geografica x (caratterizzata dalle coordinate xu e xv) all’interno dell’area S, è, secondo la definizione data in precedenza, la Variabile Regionalizzata (VR) che descrive il fenomeno in esame. Essa è chiaramente una funzione deterministica. Consideriamo ora un particolare punto di S di posizione x0. In esso si può definire una variabile aleatoria (VA) continua Z(x0), cioè una variabile che assume dei valori numerici appartenenti ad un certo intervallo secondo una legge di densità di probabilità f0 (Z) (vedi fig. 3.7).

9

Fig. 3.7 - Descrizione di una Funzione Aleatoria nel dominio S.

In queste condizioni il valore z(x0) della profondità della falda nel punto x0 può essere considerato come una realizzazione della VA Z(x0). Così come in x0, in ogni altro punto x di S può essere definita una VA Z(x). Allora l’insieme di tutte le VA definite in S costituisce una Funzione Aleatoria (FA). La FA Z(x) sarà caratterizzata dall’insieme di tutte le funzioni di distribuzione multivariabili:

FZ 1,Z 2,...,Zk ( z1, z 2,..., zk ) = prob{Z ( x1) < z1, Z ( x 2) < z 2,..., Z ( xk ) < zk ,} che si possono definire in S per qualsiasi intero k e per qualsiasi configurazione dei k punti: x1 , x2, ... xk,. Questo insieme di funzioni, data la natura spaziale del fenomeno che si vuole modellizzare, costituisce la legge spaziale della FA Z(x). Nel caso di indipendenza, a due a due, delle variabili Z(x1), Z(x2),…,Z(xk), la legge spaziale di Z(x) è costituita dall’insieme di tutte le funzioni di distribuzione monovariabili:

F Z ( x ) ( z ) = prob {Z ( x ) < z } ∀ x ∈ S Con riferimento all’esempio della falda, considerando i valori delle profondità z(x) come particolari realizzazioni delle VA definite in tutti i punti x di S, possiamo affermare che l’interpretazione probabilistica del fenomeno di studio (e cioè l’approccio probabilistico) consiste nel considerare la VR z(x) come una particolare realizzazione della FA Z(x). L’interpretazione probabilistica dei fenomeni spaziali, od il ricorso, come si usa dire, ad un modello topo-probabilista, non corrisponde ad una particolare concezione della realtà: questo modello non è l’immagine di nessuna realtà fisica, è un intermediario di calcolo, una forma di collocarsi all’interno di un quadro imperativo, in cui è possibile usare gli strumenti del calcolo delle probabilità, senza i quali, dotati semplicemente di modelli empirici o deterministici, sarebbe difficile fare molta strada verso la soluzione dei problemi. 10

E’ innegabile quindi il vantaggio dell’approccio probabilistico, ma le condizioni sono che è necessario disporre di modelli. Si pone subito quindi un problema metodologico: il problema dell’inferenza statistica, cioè della stima dei parametri del modello mediante i dati di una campionatura del fenomeno. Il modello più completo sarebbe costituito dalla legge spaziale sopra definita. Ma tale legge è impossibile da stimare a partire da una realizzazione unica, per di più conosciuta solo in un limitato numero di punti: sarebbe come pretendere di stimare con un solo risultato la legge di probabilità associata al lancio di un dado, e quindi capire, per es., se un dado è truccato o no. E’ necessario quindi fare delle ipotesi per rendere più semplice il modello. Ma è evidente che la formulazione delle ipotesi deve essere tale che il modello risultante abbia dei requisiti minimali, quali per esempio: • il modello deve poter essere stimabile nei suoi parametri, considerando il tipo e la quantità di informazioni generalmente disponibili per tale operazione; • il modello deve, come strumento, essere sufficiente per effettuare le operazioni più frequentemente richieste nello studio dei fenomeni di interesse (stime di variabili, stime di probabilità, simulazioni); • a prescindere dalle operazioni che devono essere effettuate, il modello deve poter esprimere, in forma qualitativa e quantitativa, la variabilità spaziale del fenomeno di studio, ai fini della sua comprensione. In relazione al tipo di operazioni, tra quelle sopra ricordate, che con l’ausilio del modello devono essere condotte, accenniamo fin da ora ai requisiti che esso deve possedere. Per ciò che concerne le stime (predizioni in termini di statistica classica) le più diffuse sono quelle lineari. Per la loro implementazione, sarà mostrato, è sufficiente, dell’intera legge spaziale, la conoscenza dei primi due momenti della FA; è ad essi quindi che le proprietà del modello saranno riferite. Per ciò che riguarda le simulazioni e le stime non lineari, quest’ultime comprendenti le stime di probabilità, le proprietà del modello riguarderanno non solo i primi due momenti ma anche la legge bivariabile delle coppie, cioè la legge spaziale appena introdotta limitata a k = 2.

3.3 Modelli di Base

Come è stato accennato nel paragrafo precedente il ricorso ad un modello probabilistico solleva un problema metodologico: la scelta del modello e la stima dei suoi parametri. Il problema è tutt’altro che banale, considerando che, da un lato, si è di fronte ad una vastità di modelli di Funzioni Aleatorie e, dall’altro, si dispone solo di una realizzazione unica, per di più nota solo in pochi punti del campo. Per superare le difficoltà poste da questo problema sembra opportuno porre dei vincoli ai modelli, sì da ridurre la famiglia di quelli praticamente utilizzabili.. A questo punto è utile fare la considerazione seguente:

11

la maggior parte dei problemi con cui normalmente si ha a che fare, e che sono quelli a cui si fa riferimento in queste dispense, sono problemi locali, cioè problemi, generalmente di stima, che interessano una piccola porzione del campo che comprende il supporto della stima, cioè l’entità geometrica da stimare, e l’insieme dei dati con cui effettuare la stima stessa. Questa porzione di campo, così configurata, la chiameremo d’ora in poi, con una parola presa a prestito dal francese, vicinaggio. La nozione di vicinaggio giuoca un ruolo metodologico molto importante nella Geostatistica. Poichè in generale non ci si limita ad un solo problema locale, ma si tende a fare delle stime su tutto il campo, questo vicinaggio, quasi identico per dimensioni e configurazione, è destinato a percorrere tutto il campo: diventa un vicinaggio mobile. La considerazione precedente fa sorgere una esigenza legittima: lavorando in vicinaggio mobile si vorrebbe vedere nella ripetizione spaziale legata alla traslazione del vicinaggio, un sostituto di una molteplicità di realizzazioni della FA. E ciò è tanto più leggittimo quanto più è lecito pensare ad un modello che non coinvolga tutto il campo, ma che sia definito al livello di vicinaggio. Ecco configurarsi la necessità di sollecitare per la FA condizioni di stazionarietà. Una ulteriore spinta a questo tipo di sollecitazione viene dal fatto che si dispone di una realizzazione unica del modello.

3.3.1 Modelli Stazionari Stazionarietà Strictu Sensu

La stazionarietà strictu sensu di una FA è una proprietà che si definisce molto chiaramente: è l’invarianza per traslazione della legge spaziale del processo aleatorio. Più esplicitamente: comunque siano fissati un numero intero k > 0 ed un insieme di k punti in S aventi posizione: {x1, x 2,..., xk } , i due vettori aleatori: {Z ( x1), Z ( x2),..., Z ( xk )} e {Z ( x1 + h), Z ( x2 + h),..., Z ( xk + h)} con h vettore qualsiasi , hanno la stessa legge o, più specificamente, hanno la stessa funzione di distribuzione k-variabile. La legge spaziale, anche nel quadro di un’ipotesi molto forte come la stazionarietà, continua a rimanere, per gli stessi motivi specificati nel paragrafo predente, un riferimento teorico con scarso o nullo significato pratico (se non nel caso di multigaussianità). Poichè, come è stato detto nel paragrafo precedente, la maggior parte delle operazioni geostatistiche richiedono solo la conoscenza dei primi due momenti della FA, è a questi ultimi che la proprietà della stazionarietà deve essere riferita. Prima di esaminare questa questione introduciamo i momenti del primo e del secondo ordine di una FA. Momento Primo Sia S il campo di indagine. In accordo con l’interpretazione probabilistica in ogni punto x ∈ S è definita una VA Z(x). Il suo momento primo: E Z ( x) = m ( x )

12

se esiste è generalmente funzione di x. Momento secondo Si considerino due punti x1 e x2 entrambi appartenenti ad S , rispettivamente di coordinate x1u, x1v e x2u, x2v. La covarianza tra le variabili aleatorie Z(x1) e Z(x2):

Cov( x1, x 2) = C ( x1, x 2) = E{[Z ( x1) − m( x1)]⋅ [Z ( x 2) − m( x 2)]} = E[Z ( x1) ⋅ Z ( x 2)] − m( x1) ⋅ m( x 2) è generalmente funzione delle posizioni x1 e x2. Variogramma E’ la funzione più comune della Geostatistica, usata nelle applicazioni principalmente per caratterizzare la variabilità spaziale di un fenomeno regionalizzato. Introdotto empiricamente nel paragrafo 3.2 sarà ora presentato nel quadro probabilistico. Sia S il dominio in cui è definita la FA Z(x) e siano x0 e x0+h una coppia di punti appartenenti ad S e distanti h1 (fig.3.4) La differenza tra Z(x0 ) e Z(x0+h) definisce una nuova VA detta accrescimento o incremento: [Z(x0+h) - Z(x0)]. La sua semi-varianza, è per definizione il variogramma:

γ ( x 0, h ) =

1 Var{[Z ( x 0 + h) − Z ( x 0)]} . 2

Esso in generale è funzione di h e del punto di appoggio x0 . Modelli Stazionari di Ordine 2

Un modello di FA si dice essere stazionario di ordine 22 quando sono verificate le due seguenti condizioni:

• il momento primo esiste ed è invariante rispetto alla posizione x:

1 Il simbolo di vettore indica che la distanza è considerata tenendo conto anche della direzione individuata dalla coppia di punti. 2 Nella letteratura geostatistica è frequente l’uso del termine stazionario senza l’ulteriore specificazione “di ordine 2” o “strictu sensu”. In questi casi la stazionarietà, o eventualmente la non stazionarietà, sono generalmente da intendendersi di ordine 2, rimanendo comunque inteso che il corretto significato lo si deduce sempre dal contesto.

13

E Z ( x ) = m;

(3.0)

• la covarianza C ( x 1, x 2 ) esiste e non dipende dalla particolare posizione di x1 e x2 ma dipende solo dalla loro distanza h: C ( x 1, x 2 ) = C ( x 1, x 1 + h ) = C ( h ) .

(3.1)

L’esistenza della funzione covarianza implica che per x 2 → x1 , cioè per h → 0 il valore della covarianza C(0) esiste ed ha valore finito. Esso corrisponde alla varianza di Z(x), detta anche varianza a priori, che è quindi anch’essa invariante per traslazione: C ( x 1, x 1 ) = Var ( x 1) = C ( 0 ) .

Quando sono verificate le due condizioni precedenti la funzione C(h) è chiamata funzione covarianza. Essa esprime la correlazione tra le variabili Z(x) e Z(x+h) in funzione della distanza h dei loro punti di appoggio: C ( h ) = E Z ( x + h ). Z ( x ) − m2 Proprietà della covarianza:

• è una funzione pari: C(h) = C ( −h) ;

• la covarianza nel punto h = 0 assume sempre valori positivi: C ( 0 ) = var Z ( x ) ≥ 0 ;

• vale la disuguaglianza di Schwartz: C ( h ) ≤ C ( 0)

( 3.2 )

Poichè, come è ragionevole pensare, la correlazione tra le variabili Z(x) e Z(x+h) tende ad indebolirsi con l’aumento della distanza tra i due punti di appoggio, si ha che la funzione C(h) tende a decrescere con |h| = h, fino anche a potersi annullare se le due variabili diventano indipendenti.

14

a: range

C(h) C(0)

a

h

Fig. 3.8. Andamento generale della funzione covarianza.

Se ciò ha luogo la distanza alla quale le due variabili diventano indipendenti si chiama portata (range nella terminologia inglese e portée in quella francese) (fig.3.8) Esaminiamo ora il comportamento della funzione variogramma nel quadro di un modello stazionario di ordine 2. Partiamo dalla definizione:

1 2

γ ( x , h ) = Var Z ( x + h ) − Z ( x ) =

{

( 3.3)

}

1 2 2 E [Z ( x + h) − Z ( x)] − {E [Z ( x + h) − Z ( x)]} 2

Poichè, per ipotesi, il momento primo della FA Z(x) è invariante per traslazione, si ha che: E Z ( x + h) − Z ( x) = 0.

Quindi la ( 3.3 ) diventa:

γ (x , h ) =

1 E 2

{[ Z

( x + h ) − Z ( x )]

2

}

.

(3.4)

Il variogramma alla distanza h coincide, a meno del fattore ½ con la media degli accrescimenti quadrati di Z(x) di entità h. Ancora la ( 3.4 ) sviluppata secondo la relazione ( 1.20) diventa:

γ ( x, h ) =

1 {Var[Z ( x + h)] + Var[Z ( x)] − 2Cov[Z ( x + h), Z ( x)]} 2

(3.5)

Ora, sempre in virtù delle ipotesi di stazionarietà di ordine 2, la varianza e la funzione covarianza sono invarianti per traslazione e quindi la ( 3.5 ) diventa:

15

γ (h ) =

1 { 2Var [Z (x )] − 2C (h )} , 2

da cui, riscrivendo più chiaramente, si ottiene l’importante relazione:

γ ( h ) = C ( 0) − C ( h )

( 3.6 )

Questa relazione mostra che la funzione variogramma è strettamente legata alla funzione covarianza e pertanto si può affermare che: • anche γ(h) è invariante per traslazione; • è indifferente utilizzare, come strumento di lavoro, la covarianza o il variogramma.

C (0 )

g (h)

C (h)

R ange

h

Fig. 3.9- Funzione covarianza e funzione variogramma.

Dalla ( 3.6 ), tenendo conto della relazione di Schwartz ( 3.2 ), si deduce che la funzione γ(h) è necessariamente limitata; il limite superiore è rappresentato dalla varianza a priori (fig. 3.9).

3.3.2 Modelli non Stazionari Le FA non stazionarie sono quelle che non rispettano le condizioni (3.0) e (3.1). Più esplicitamente sono quelle che soddisfano anche ad una delle due condizioni: • la media E Z ( x ) = m( x ) non è costante nel campo; • la funzione covarianza o non esiste (perchè non esiste la varianza) o non è invariante per traslazione. Prima di esaminare le forme con cui le FA non stazionarie vengono affrontate e caratterizzate al fine di rendere possibile la loro manipolazione, ci soffermiamo su due tipi

16

di modelli, non stazionari secondo la definizione data, che in maniera molto immediata si riconducono allo stesso formalismo delle FA stazionarie. Modelli quasi Stazionari Una FA è detta essere quasi-stazionaria quando la sua media E[Z(x)] = m(x), pur non essendo costante su tutto il campo, vi varia molto debolmente, sì da poter essere considerata costante all’interno di domini di dimensioni non inferiori a quelle del vicinaggio di lavoro. Vediamo come si esprime un modello quasi-stazionario. Consideriamo, per svincolarci dalla notazione vettoriale, di operare su uno spazio a una dimensione, o su un allineamento (uno qualsiasi) di un campo bidimensionale. Siano x0 e x0 +h due punti sull’allineamento. Il variogramma di Z(x) per la coppia Z(x0) e Z(x0+h), in base alla ( 3.3) è dato da:

{

2γ (x 0 , h ) = E [ Z (x 0 + h ) − Z (x 0 )]

2

} − [ m (x

0

+ h ) − m (x 0 )]

2

E’ facile verificare che il variogramma sopra definito è anche il variogramma di: Y(x) = Z(x) - m(x), che è una variabile, detta residuo, di media nulla. Supponiamo che tale residuo sia stazionario di ordine 2, il che equivale a supporre, essendo il residuo di media costante, che γ ( x 0 , h ) dipenda solo da h e sia limitato superiormente. Supponiamo ancora che m(x) sia una funzione, per semplicità, di tipo lineare: m( x ) = a 0 + a 1 x . l’espressione :

{

}

1 2 E [Z ( x 0 + h) − Z ( x 0)] , 2 che coinciderebbe con il variogramma di Z(x) se E Z ( x ) fosse costante, risulta essere: 1 2 E [ Z ( x 0 + h ) − Z (x 0 )] = γ (h ) + a 12 h 2 ( 3.7 ) 2

{

}

Ora, se il coefficiente a1 è molto piccolo, i termini di destra della espressione precedente, dove il primo rappresenta un variogramma limitato (fig. 3.10 a), il secondo una funzione parabolica (fig. 3.10 b), sommati danno luogo all’andamento di fig. 3.10 c, che può essere considerato come l’andamento di un variogramma di una FA Z(x) stazionaria entro i limiti fissati dalla distanza b alla quale la curva comincia ad impennarsi. Tale distanza definisce il dominio di stazionarietà della FA Z(x). In altri termini Z(x) può essere considerata stazionaria di ordine 2 per vicinaggi di dimensioni ≤ b ed il suo variogramma è espresso da:

γ (h) ≅

{

1 2 E [Z ( x 0 + h) − Z ( x 0)] 2

} 17

In pratica si è considerato, sempre per h ≤ b, a 1 2 h 2 ≅ 0 . E’ importante notare che la illustrazione del concetto di quasi-stazionarietà è stata effettuata servendosi della funzione variogramma; non sarebbe stato possibile usando la covarianza. E’ questa una delle ragioni per cui nella geostatistica applicata si opera in principio con la funzione variogramma.

g (h)

a)

h 2

m (h)

b)

h

c)

g *(h)

h

Fig. 3.10 - Andamento del variogramma sperimentale per modelli quasi-stazionari.

Modelli intrinseci

18

Progressiva

Può succedere che una FA Z(x) non stazionaria sia tale che i suoi incrementi lo siano. Una situazione del genere è classica e corrisponde al moto browniano, vale a dire al moto disordinato di una particella provocato dagli urti con una grande quantità di altre particelle in movimento costituenti un liquido. Il risultato di una simulazione del moto browniano a una dimensione, su 6000 unità temporali, è riportata in figura 3.11 e consiste nell’andamento, in funzione del tempo, della progressiva che identifica la posizione della particella.

80 60 40 20 0 -20 -40 -60 -80 0

1000

2000

3000

4000

5000

6000 Tempo

Fig. 3.11 – Andamento delle progressive (simulate) relative al moto (browniano monodimensionale) casuale di una particella a seguito degli innumerevoli urti subiti ad opera di altre particelle.

Come si evince dal grafico la variabile presenta un andamento senza regole che non consente di fare alcuna considerazione su media e varianza. Gli incrementi invece sono stazionari, nel senso che essi fluttuano attorno un valore costante (zero in particolare) e con fluttuazioni di ampiezza paragonabile. Nella figura 3.12 si riporta l’andamento di 500 incrementi a passo temporale 3 e 6 calcolati nella parte iniziale, nella parte di mezzo e nella parte finale della simulazione. L’andamento di fig. 3.11 è modellizzabile con una FA intrinseca.

19

0

100

200

300

400

Incrementi di progressiva (Dt=6)

Incrementi di progressiva (Dt=3)

(1)

4 3 2 1 0 -1 -2 -3 -4

8 6 4 2 0 -2 -4 -6 -8 0

500

100

200

100

200

300

400

Incrementi di progressiva (Dt=6)

Incrementi di progressiva (Dt=3)

(2)

4 3 2 1 0 -1 -2 -3 -4 0

500

0

300

400

500

Incrementi di progressiva (Dt=6)

Incrementi di progressiva (Dt=3)

(3)

200

500

100

200

300

400

500

Tempo

4 3 2 1 0 -1 -2 -3 -4 100

400

8 6 4 2 0 -2 -4 -6 -8

Tempo

0

300

Tempo

Tempo

8 6 4 2 0 -2 -4 -6 -8 0

Tempo

100

200

300

400

500

Tempo

Fig. 3.12 – Andamento delle progressive (simulate) a passo temporale 3 (parte sinistra) e 6 (parte destra) calcolati rispettivamente sui valori 1-500 (1), 2500-3000 (2) e 5500-6000 (3) della simulazione.

Una FA si dice intrinseca in S (e si indica FAI) se per ogni vettore h i corrispondenti accrescimenti sono stazionari di ordine 2, vale a dire: dati due punti di posizione x0 e x0 + h l’accrescimento: Z(x0 ) - Z(x0 + h) ammette momenti ordine uno e due di valore finito ed invarianti per traslazione, cioè non dipendenti dal punto di appoggio x0 ma solo dalla posizione h. Le due condizioni si esprimono come segue: segue: = [m(h) E Z ( x 0 + h ) − Z ( x 0 ) ] = m(h)

( 3.8 )

D 2 Z ( x 0 + h ) − Z ( x 0) = 2 γ ( h )

( 3.9 )

m(h) rappresenta la deriva di Z(x) che è necessariamente lineare. Infatti la (3.8) si può anche scrivere:

[

]

’

E Z ( x 0 + h + h ') − Z ( x 0 ) = m( h + h )

e anche

20

(3.10)

[

]

’ E Z ( x 0 + h + h ') − Z ( x 0 + h ) = m(h )

(3.11)

Sottraendo la (3.8) dalla (3.10) e tenendo conto della (3.11) si deduce la linearità di m(h): m(h+h’) - m(h) = m(h’). Le FAI cui si farà riferimento nel seguito saranno quelle senza deriva e cioè con m(h) = 0. Questa restrizione, che nell’ambito della Geostatistica di base è utile per la esemplificazione delle procedure, deve essere verificata. In un ambito più generale, quello in cui si studiano i fenomeni non stazionari, tale restrizione è invece fondata. Si precisa che una FAI senza deriva non significa una FA di media costante, perché, per definizione, una FAI non fa nessun riferimento alle VA Z(x) ma solo agli incrementi. L’andamento di fig. 3.11 non evidenzia una fluttazione attorno ad un valore costante, eppure gli incrementi sono di media nulla. Nella (3.9) γ(h) rappresenta la funzione variogramma già introdotta all’inizio del paragrafo. Il variogramma, appunto perchè espressione di una FA intrinseca, è anche chiamata nella terminologia geostatistica funzione intrinseca. Un modello intrinseco è chiaramente meno esigente di un modello stazionario: le condizioni che lo definiscono sono meno restrittive. Si osserva che un modello stazionario è necessariamente intrinseco: infatti nel paragrafo 3.3.1 è stato dimostrato che se un modello è stazionario esso ammette variogramma, che risulta essere legato alla funzione covarianza dalla relazione (3.6). Si ribadisce che un modello intrinseco non ammette obbligatoriamente funzione covarianza: si ricorda anzi che il modello intrinseco è stato proprio introdotto per trattare FA prive di covarianza. Un modello (o una FA) intrinseco ma non stazionario è chiamato strettamente intrinseco. Senza entrare nel merito di cose che saranno esaminate più avanti, è importante tener presente in questa fase che nella stima lineare il formalismo di calcolo per le FA intrinseche è identico a quello delle FA stazionarie, a patto di sostituire la funzione covarianza con la funzione variogramma cambiata di segno.

I due modelli sopra presentati, molto ricorrenti nella pratica geostatistica, soprattutto il primo, sono modelli non stazionari molto particolari, che, come si è visto, sono stati facilmente ricondotti a situazioni di stazionarietà: il primo restringendo la stazionarietà su aree ristrette, il secondo cercando la stazionarietà dell’incremento. Anche per le altre situazioni di non stazionarietà, quelle più proprie, ma anch’esse molto frequentemente ricorrenti nella modellizzazione geostatistica, si seguirà un cammino analogo: attraverso il ricorso a particolari forme ci si ricondurrà sempre a un qualcosa di stazionario. Due sono i tipi di modelli proposti nella geostatistica non stazionaria: modello con deriva e modello intrinseco di ordine k: sono due facce della stessa medaglia (non stazionaria); entrambi filtrano dal trattamento probabilistico tutto ciò che non è caratteristica intrinseca del fenomeno; il secondo modello è una generalizzazione del modello intrinseco.

21

Senza entrare nei particolari, diamo qui di seguito un breve cenno sugli approcci per completare il quadro sui modelli delle Funzioni Aleatorie.

Modello con Deriva Questo tipo di modello trae spunto dall’osservazione di alcuni fenomeni naturali: quelli che mostrano una componente di variabilità a scala regionale o una variabilità che ingloba un trend, vale a dire una variazione sistematica della variabile più o meno accentuata, in un senso o nell’altro. Empiricamente questa tendenza potrebbe essere descritta dall’andamento della media (della variabile di interesse) su vicinaggi mobili all’interno dell’area di indagine. Il punto di vista del vicinaggio determina la scala alla quale bisogna osservare il fenomeno ai fini dello studio geostatistico, perchè, è evidente, certe frequenze appaiono o non appaiono a seconda della scala di osservazione. Se un fenomeno quindi manifesta con evidenza, alla corretta scala di osservazione, un progressivo aumento o diminuzione, vuol dire che bisogna scegliere una FA modello con media variabile: E Z ( x ) = m( x ) ; quest’ultima nel linguaggio geostatistico è chiamata deriva (dérive nella terminologia francese e trend in quella inglese), termine già introdotto per caratterizzare una FAI. Pertanto nel modello con deriva è supposto che la FA oggetto di studio Z(x) possa essere considerata, in ogni punto x, come la sovrapposizione di due componenti: • la deriva m(x), che costituisce la componente deterministica; • il residuo Y ( x ) = Z ( x ) − m ( x ) , che costituisce la parte aleatoria. Si ha ovviamente che: Z ( x ) = Y ( x ) + m( x ) . Questa dicotomia, decomponendo la variabilità totale in due componenti, una deterministica e l’altra stocastica, con un criterio che, come è stato detto sopra, è influenzato dalla scala di osservazione del fenomeno, sembra introdurre un forte elemento di soggettività. In verità, spesso questa dicotomia risponde ad una esigenza naturale, e ciò accade quando la deriva è espressione di una componente regionale, regolare, di una tendenza che ha un significato fisico rispetto al fenomeno studiato. Generalmente la deriva viene modellizzata con una funzione polinomiale: m( x ) =

al f l ( x ) , l

dove al sono i coefficienti del polinomio e fl(x) monomi di grado l. Il residuo Y(x), essendo esso a media costante (perchè nulla), viene modellizzato con una FA o stazionaria o intrinseca a seconda che ammetta covarianza o solo variogramma. In un caso avremo un modello con deriva e residuo stazionario, nell’altro un modello con deriva e residuo intrinseco. Al variogramma e la covarianza del residuo, espressi nella forma consueta: 22

C ( h) = E Y ( x ) ⋅ Y ( x + h ) ;

{

2γ (h ) = E [Y (x + h ) − Y (x )]

2

}

è attribuito l’appellativo di sottogiacenti, perché la loro inferenza, quando è possibile, non è diretta.

Modello Intrinseco di ordine k La definizione di questo modello è il risultato della generalizzazione del procedimento seguito per la definizione del modello strettamente intrinseco: allora si è operato considerando in ogni punto x invece della variabile diretta Z(x) l’accrescimento in x relativo ad una distanza orientata h: Z(x+h) - Z(x). In questo modello la nozione di accrescimento è generalizzata ad una combinazione lineare di ordine k, cioè ad una combinazione lineare: n

Z Λk ( x ) =

λiZ ( x i )

(3.12)

i =1

che, affinché possa essere stazionaria, deve filtrare dei polinomi di grado k, cioè deve soddisfare alle seguenti condizioni: n

λif l ( x i ) = 0

(l = 0, 1,...,k)

(3.13)

i =1

Nella (3.12) {x1, x 2,..., xn} sono le posizioni di n punti ∈S ; {Z ( x1), Z ( x 2),..., Z ( xn)} sono le VA ad essi associate; Λk un insieme di numeri reali {λ 1, λ 2,..., λn} , che costituiscono i coefficienti della combinazione lineare.

Nella (3.13) f l ( x i ) sono i monomi di grado l delle coordinate degli n punti. Nel caso che il campo S sia bidimensionale, indicando con {xu1, xu 2,..., xun} e {xv1, xv 2,..., xvn} le coordinate degli n punti, per un ordine per es uguale a 2, la (3.13) equivale ai seguenti tre gruppi di condizioni: I° gruppo:

II°gruppo:

{

i

λi = 0

λixui = 0 λixvi = 0 i

i

λix ui = 0 λix vi = 0 i λixuixvi =0 2

III°gruppo:

i

2

i

23

Per un ordine k = 1 le condizioni sono quelle del primo e del secondo gruppo, mentre per un ordine k = 0 le condizioni sono solo quelle del primo gruppo. Possiamo ora enunciare le proprietà del modello. Una intrinseca di ordine k se la FA : Z Λk ( x ) =

n

FA

Z(x) si dice essere

λiZ ( x i ) è stazionaria di ordine 2. In tal

i =1

caso la combinazione Z Λk ( x ) è detta combinazione lineare autorizzata di ordine k. A questo punto non è difficile riconoscere nell’accrescimento stazionario Z(x+h) Z(x) una combinazione lineare autorizzata di ordine 0. Infatti in questo caso i coefficienti dei termini della combinazione lineare sono 1. e -1. e rispettano la condizione del primo gruppo.

3.3.3 Considerazioni sulla Scelta dei Modelli

Dopo questa disamina sui modelli di FA che ricorrono nelle applicazioni della Geostatistica inevitabilmente si pone il problema della scelta del modello con interrogativi del tipo: con quale criterio, con quali dati, con quali tecniche, con che conseguenze. A tutte queste domande verrà risposto specificamente nel successivo paragrafo 3.4 per ciò che concerne una parte dei modelli, e nel Cap. 4 per la restante parte, atteso che tutto lo sviluppo di queste dispense contribuisce in maniera generale a chiarire questo problema. In questo punto ci limitiamo a qualche considerazione di principio. Innanzi tutto va chiarito il rapporto fenomeno/modello. Le proprietà che definiscono un modello sono delle proprietà matematiche e sono da attribuire unicamente alla FA. Non si possono attribuire delle proprietà matematiche ai fenomeni. Con riferimento a una delle scelte più ricorrenti nella pratica geostatistica, se adottare un modello stazionario o un modello non stazionario per lo studio di un fenomeno, possiamo efficacemente citare P. Chauvet (“Aide mémoire de Géostatistique linéaire”, pag.89) “... una variabile regionalizzata (quindi un fenomeno) non è, in sé, nè stazionaria, nè non stazionaria. Per essere precisi bisognerebbe dire tale variabile regionalizzata, per un dominio dato e per una data scala di lavoro, può ragionevolmente essere (o non essere) modellizzata da una Funzione Aleatoria stazionaria. E si sottolinea ancora una volta che l’avverbio ragionevolmente, che può scandalizzare gli statistici, potrebbe essere presa come la pietra angolare del procedimento del geostatistico. Così nel caso preciso che ci interessa, non c’è vera risposta alla questione stazionario o non stazionario, ma piuttosto una assunzione di responsabilità da parte del geostatistico, fatta con cognizione di causa, nel migliore accordo possibile coi dati, seguendo sempre il principio dell’economia, e cosciente che comunque alla fine ci sarà sempre un rischio”. Bisogna precisare che in geostatistica non si ha alcuna avversione pregiudiziale per i test, ma si deve tener presente che questi test necessitano di disporre già di un modello interamente specificato. Anzi, proporre dei test di stazionarietà esigerebbe, a livello teorico, una specificazione del modello probabilista molto più completa di quanto non sia richiesto ulteriormente dalla Geostatistica Lineare e di quanto realisticamente si possa sperare di raggiungere.

24

Dal punto di vista pratico la scelta del modello non presenta poi particolari difficoltà se si ha un pizzico di esperienza e si dispone di software geostatistici adeguati. Il più delle volte il modello da adottare viene suggerito in maniera inequivocabile dal comportamento dei dati rispetto a certe operazioni. E ciò non è casuale dato che i modelli proposti traggono spunto dalla realtà. Quando poi sorgono delle difficoltà ricordiamo che esiste sempre un criterio empirico che aiuta a fare scelte efficaci: quello di scegliere il modello sulla base dei risultati di una validazione incrociata.

3.4 Analisi ed interpretazione dei variogrammi sperimentali

Con riferimento ai modelli di FA esaminati nel paragrafo precedente, si ha che il modello stazionario, il modello quasi-stazionario ed il modello intrinseco possono essere tutti e tre descritti nella stessa forma: tramite la funzione variogramma. In tutti e tre i casi infatti, per le proprietà del modello, la funzione variogramma esiste ed è invariante per traslazione, con una restrizione per il modello quasi-stazionario che limita queste proprietà al dominio di stazionarietà. Per questi modelli si ricorda che il variogramma ha la seguente forma:

γ (h) =

{

1 2 E [Z ( x + h) − Z ( x)] 2

}

(3. 14)

Questa espressione si presta ad una stima diretta a partire dai valori misurati e ciò rende piuttosto semplice l’ inferenza della funzione variogramma. Questo paragrafo sarà dedicato alla stima del variogramma e della sua interpretazione nell’ambito dei modelli citati. Per gli altri modelli, che sono quelli tipicamente non stazionari, la definizione ed il riconoscimento del modello richiedono procedure diverse che saranno descritte nel cap. 4, dedicato totalmente ai fenomeni “non stazionari”. Quando si affronta uno studio geostatistico, la prima fase di esso, dopo una necessaria analisi preliminare consistente nella elaborazione delle statistiche elementari e finalizzata a prendere conoscenza dei dati ed alla loro verifica, è dedicata al calcolo ed alla interpretazione dei variogrammi sperimentali. Uno dei risultati dell’interpretazione è la identificazione del modello di FA: se stazionario, quasi-stazionario, intrinseco, o qualcosa altro di non stazionario. Spesso l’ identificazione del modello, come sarà mostrato più avanti, avviene in maniera molto chiara e inequivocabile. Con riferimento ai modelli di FA introdotti ed ai metodi che saranno descritti in queste dispense si ricorda che, ad eccezione del cap. 6 dedicato al trattamento dei fenomeni non stazionari, il resto dei capitoli tratterà le metodologie per modelli stazionari e quasi-stazionari. Tutto quanto è svolto nell’attuale cap. 3 vale anche per i modelli intrinseci.

3.4.1 Calcolo dei variogrammi sperimentali

25

La stima della funzione variogramma viene effettuata sulla base dei dati provenienti dal campionamento del fenomeno oggetto di studio. Si consideri per esempio la fig. 3.13 in cui è schematizzato un campo di lavoro ed i punti di prelievo, disposti a maglia regolare, di campioni di suolo su cui è stato misurato il contenuto in metalli pesanti, tra essi il Cadmio (Cd). Si vuole effettuare il calcolo del variogramma sperimentale. Si tratta di stimare l’espressione (3.14), che per FA stazionarie, quasi-stazionarie e intrinseche, rappresenta appunto il variogramma. Il suo stimatore, come ogni altro stimatore considerato in queste dispense, viene convenzionalmente notato contrassegnando con un asterisco l’entità da stimare:

{

}

1 γ ∗ (h) = Stima E [Z ( x + h) − Z ( x)]2 . 2

(3. 15)

Fig. 3.13 - Disposizione di un piano di campionatura a maglia regolare.

Per le FA sopra menzionate, data la stazionarietà dell’ incremento Z ( x + h) − Z ( x ) , in assenza di deriva γ ∗ ( h ) è dato dalla media aritmetica degli incrementi quadrati misurati a distanza h. Siano infatti r il lato della maglia e xi , xi + r la posizione di due punti distanti r e allineati secondo la direzione X, il variogramma sperimentale alla distanza r nella direzione X è dato da:

γ ∗ (r ) =

Nr i =1

z ( xi + r ) − z ( xi )

2

dove N r è il numero di coppie di punti aventi tra di loro distanza r e allineati secondo la direzione X. Alla stessa maniera si può calcolare, secondo la stessa direzione X il variogramma sperimentale per h=2r, 3r; analogamente il tutto è ripetibile per la direzione Y. Per la direzione a 45° è possibile calcolare γ ∗ ( h ) per i valori di h: r 2 ,2r 2 ,3r 2 . . Non così immediato è il calcolo dei variogrammi sperimentali quando la maglia è irregolare, come per esempio lo schema di campionature riportato in fig. 3.14 a. Supponiamo di voler calcolare il valore del variogramma ad una distanza r lungo la direzione φ (fig. 3.14 b). Il contributo al calcolo del variogramma suddetto deriva da tutte le coppie di campioni che si possono individuare distanti r e allineate secondo la 26

direzione φ. Può accadere che non vi sia alcuna coppia così disposta. E’ allora necessario introdurre una tolleranza angolare ∆φ sulla direzione ed una tolleranza ∆r sulla distanza; così facendo tutte le coppie di campioni aventi distanza compresa tra r∆r e r+∆r e allineate secondo una direzione conpresa tra φ-∆φ e φ+∆φ contribuiscono al calcolo del variogramma r ,φ. Dal punto di vista operativo è immediato verificare, con semplici considerazioni geometriche, se la coppia di campioni individuata dai punti p1 e p2 di coordinate rispettivamente (u1 , v1 ) e (u2 , v2 ) contribuisce al calcolo del variogramma per la direzione φ ± ∆φ e la distanza r ± ∆r.

Fig. 3.14 - Disposizione di un piano di campionatura a maglia irregolare.

Infatti , indicando con d la distanza tra i due punti: d = sqrt ((u2 -u1)∗∗2+(v2 -v1)∗∗2 ), con con C1 = d / (u2 -u1) e C2 = d /(v2 -v1) i coseni direttori della direzione p1 - p2 e con D1 = cos(φ) e D2 = sin(φ) i coseni direttori della direzione φ, il coseno dell’angolo α , formato dalla congiungente i punti p1 , p2 e dalla direzione φ, è dato da: cos(α ) = C1∗D1 + C2∗D2 . Quindi la coppia contribuirà al calcolo del variogramma se: cos(α ) ≥ cos(∆φ ) e r-∆r ≤ d ≤ r+∆r. I valori delle tolleranze ∆φ e ∆r da adottare dipendono ovviamente dalla quantità di campioni di cui si dispone: più essi sono di numero elevato e più piccole possono essere le tolleranze, consentendo con ciò maggiore precisione nel calcolo dei variogrammi. E’ comunque di uso abbastanza frequente calcolare i variogrammi sperimentali per distanze multiple di una distanza di base, chiamata passo, con un valore di ∆r pari proprio a metà del passo e secondo quattro direzioni, a novanta gradi tra di loro, con un valore di ∆φ di 22.50 ° . Questo consente di utilizzare tutte le coppie individuate dall’insieme dei campioni disponibili. In conclusione, sia partendo da una maglia regolare che da una maglia non regolare, utilizzando le tecniche di calcolo sopra illustrate, è possibile stimare per diverse distanze h il variogramma sperimentale γ∗(h), i cui valori possono essere graficati come in fig. 3.15. L’andamento di γ∗(h) in funzione di h esprime quasi sempre situazioni facilmente riconoscibili, che riguardano sia il modello di FA sia lo stile della variabilità spaziale del fenomeno considerato. Esamineremo di seguito questi aspetti 27

facendo riferimento a casi di studio e fornendo di volta in volta i presupposti teorici dell’interpretazione.

Fig. 3.15 - Rappresentazione di un variogramma sperimentale.

3.4.2 Identificazione del modello di FA

Quando il variogramma sperimentale, sia pur con delle fluttuazioni, si attesta su di un valore che poi rimane costante, e tale valore è pressappoco coincidente con la varianza empirica, allora questo è un segno che il fenomeno in esame può essere descritto da un modello stazionario (di ordine 2). Si ricorda che per tale modello il momento primo ed il momento secondo sono invarianti per traslazione e che la funzione variogramma è legata alla funzione covarianza dalla relazione seguente:

γ ( h ) = C ( 0) − C ( h ) C(0) è il valore di soglia che corrisponde alla varianza della FA Z(x). Nelle tre figure seguenti riportiamo altrettanti variogrammi sperimentali, presi dalla letteratura, in cui appare molto evidente l’esistenza di un valore di soglia che rimane poi costante. In ogni figura accanto al variogramma sperimentale è tracciato il variogramma modello, cioè la funzione analitica aggiustata sui punti sperimentali; l’aggiustamento del variogramma sarà trattato nel paragrafo 3.5. La fig. 3.16 mostra il variogramma secondo la direzione verticale, costruito sui dati di 40 perforazioni, della variabile saturazione in olio della formazione sabbiosa canadese Athabasca. Il variogramma presenta delle oscillazioni attorno ad un valore stimato essere attorno a 19, raggiunto ad una distanza di 36 m. Nella figura per alcuni punti sperimentali è indicato il numero di coppie con cui il variogramma è stato calcolato.

28

Fig. 3.16 - Variogramma sperimentale della saturazione in olio di un giacimento petrolifero. (figura tratta da Mining Geostatistics di Journel & Huijsbregts, pag 237).

La fig. 3.17 si riferisce al variogramma sperimentale del contenuto in Cu del giacimento cileno di Los Bronces costruito sulla base di circa 4000 campioni. Il valore di soglia di 0.73 è stato considerato essere raggiunto asintoticamente; stando alla definizione, si direbbe che il variogramma ha un range infinito. La fig. 3.18 rappresenta il variogramma della variabile conducibilità elettrica, costruito sui dati trasformati (trasformazione gaussiana), misurata per il controllo della salinità su

Fig. 3.17 - Variogramma sperimentale della variabile Cu in un giacimento del Cile (da Mining Geostatistics di Journel & Huijsbregts, pag.240)

Fig. 3.18 - Variogramma della conducibilità elettrica ca in un suolo della valle del Giordano in Israele. (da ”Disjunctive Kriging in Agriculture”, di R. Webster and Oliver”, M.Armstrong(ed.), Geostatistics, vol. I)

campioni di suolo della regione di Bet Shean nella valle del Giordano in Israele. Come nel caso precedente il valore di soglia 1 è raggiunto asintoticamente.

29

Le situazioni sopra riportate sono inequivocabili, esse ricorrono abbastanza frequentemente nel trattamento dei dati spaziali: l’adozione di un modello stazionario appare quanto mai appropriata. Un’altra situazione, anch’ essa frequente nella pratica delle applicazioni, è quella in cui il variogramma sperimentale raggiunge un valore di soglia che si mantiene costante per un breve tratto e poi si impenna crescendo indefinitamente. Questo comportamento lo si può osservare nella fig. 3.19, che rappresenta il variogramma verticale di una mineralizzazione di Cu, in cui il contenuto di Cu diminuisce debolmente e sistematicamente con la profondità. L’andamento di un siffatto variogramma sperimentale indica inequivocabilmente, sulla base di quanto detto nel punto 3.3.2, un comportamento della VR modellizzabile con una FA quasi-stazionaria. Il valore di h al quale il variogramma sperimentale comincia ad impennarsi è di 100 ft, e corrisponde alla distanza al di sotto della quale il trend di variazione del contenuto in Cu è trascurabile. Operando al di sotto dei 100 ft si ha il diritto di modellizzare la variabile di studio con una FA stazionaria, il cui variogramma è rappresentato appunto dal tratto di curva al di sotto dei 100 ft. E’ importante notare che in questi casi il valore di soglia del variogramma non corrisponde, al contrario di quanto accade nel caso stazionario, alla varianza empirica dei campioni, in quanto questa, essendo calcolata su tutti i campioni, ingloba anche la variabilità dovuta alla variazione sistematica.

Fig. 3.19 - Variogramma sperimentale indicante una situazione di quasistazionarietà (fig. tratta da tratta da Mining Geostatistics di Journel & Huijsbregts, pag 44).

Facendo ancora riferimento alle considerazioni sulla quasi-stazionarietà del già citato paragrafo 3.3.2, si ha che, se la variazione sistematica della media non è più debole ma accentuata, la somma delle due componenti, il variogramma e la parabola, produce una curva che si impenna rapidamente non evidenziando alcun dominio di stazionarietà, come risulta dalla fig. 3.20: ci si trova di fronte ad una situazione non stazionaria.

30

Essendo i fenomeni naturali strutturati e spesso regolati da fattori che agiscono su direzioni preferenziali, può accadere che un fenomeno sia stazionario in una direzione e non stazionario in un’altra direzione. Si osservi, ad esempio, la fig. 3.21. In essa sono

g*(h)

2

m (h)

g(h)

h

h

Fig. 3.20 - andamento del variogramma sperimentale in presenza di un trend accentuato

riportati i variogrammi sperimentali del geopotenziale1 a 500 mbar, uno in direzione EW (longitudine), l’altro in direzione NS (latitudine). Il variogramma EW presenta un valore di soglia raggiunto ad una distanza di circa 3000 Km, denunciando una situazione di stazionarietà o meglio di quasi-stazionarietà. Il variogramma NS mostra invece un marcato andamento parabolico, che, per quanto detto sopra, è il segnale della presenza di una deriva, cioè di una sensibile variazione sistematica del geopotenziale.

Fig. 3.21 - variogrammi del geopotenziale (da Mining Geostatistics di Journel & Huijsbregts, pag. 242).

1

Il geopotenziale è una grandezza usata in meteorologia e corrisponde all’altezza H(x) in un punto x di una determinata isobara superficiale.

31

A conferma di questa interpretazione si osservino nella figura seguente gli andamenti misurati del geopotenziale lungo le due direzioni precedenti. Mentre l’andamento nella direzione EW (fig. 3.22b) mostra delle fluttuazioni attorno ad una

Fig. 3.22 - Andamento del geopotenziale nelle direzioni N-S, E-W (da Mining Geostatistics di pag. 242243).

media pressocchè costante, l’andamento nella direzione NS (fig. 3.21a) mostra un progressivo aumento del geopotenziale dal polo all’equatore. Diamo ancora un altro esempio di come passando da una direzione all’altra si passa da una situazione stazionaria ad una non stazionaria. Nella fig. 3.23 è raffigurata un’area di 100 mila ettari lungo il fiume Batanghary (ovest di Sumatra). L’area è stata campionata mediante il prelievo di circa 150 campioni, su cui sono state misurate le caratteristiche del suolo, tra cui la tessitura, che, in termine di scienza dei suoli, consiste nelle percentuali di sabbia, limo, argilla. La tessitura come è noto è legata alla morfologia dell’area. Sui dati di campionatura è stato condotto uno studio geostatistico finalizzato allo studio della variabilità del suolo. Nella fig. 3.24 si riportano i variogrammi del contenuto in sabbia calcolati su quattro direzioni. Come si può osservare il variogramma lungo la direzione SE-NW (che pressappoco corrisponde a quella dell’asse del fiume) raggiunge un valore di soglia di circa 150 ad una distanza di 15-20 Km. Nella direzione ortogonale alla precedente NE-SW il variogramma mostra invece un comportamento parabolico, che è il segno della presenza di un trend nella variabile. Nelle direzioni intermedie i variogrammi hanno un comportamento intermedio, anche se non simmetrico. In particolare si nota che la direzione E-W presenta una situazione di quasi-stazionarietà con dominio di stazionarietà di circa 5 Km, somigliando un po'alla dirzione SE-NW, mentre il variogramma nella direzione N-S somiglia più a quello della direzione NESW. Ciò si spiega molto semplicemente: l’asse del fiume non coincide esattamente con la direzione NW-SE e quindi le direzioni N-S, E-W non sono simmetriche rispetto all’asse del fiume. Accade che l’asse del fiume forma con la direzione E-W un angolo minore di quello formato con la direzione N-S.

32

Fig. 3.23 - Area di Sitiung (Est Sumatra). Principali unità geomorfiche e localizzazione dei campioni (da “soil variability of soil properties”, di G. Uehara et al., Proceedings of the workshop on soil spatial variability of the ISSS and the SSSA, Las Vegas 1984).

I due esempi precedenti non sono dei casi isolati; sono molto frequenti i fenomeni in cui una variabile di studio manifesta una tendenza sistematica ad aumentare o a diminuire, e ciò quasi sempre avviene secondo una determinata direzione. Il modello di FA da adottare è ovviamente un modello complessivo non stazionario; come e con quale approccio sarà illustrato nel par. 4.

Fig. 3.24 - Variogrammi della percentuale di sabbia calcolati su quattro direzioni del piano (elaborati da “soil variability of soil properties”, di G. Uehara et al., Proceedings of the workshop on soil spatial variability of the ISSS and the SSSA, Las Vegas 1984).

33

Non necessariamente quando il variogramma sperimentale non raggiunge un valore di soglia, esso deve avere un comportamento parabolico. E’ frequente il caso in cui l’andamento è lineare, come il variogramma mostrato nella fig. 3.25b. Esso si riferisce al variogramma dei mm d’acqua misurati dopo un evento nel bacino di Kadjemeur nel Tchad (Africa Centrale). La sottostante figura 3.25a mostra la localizzazione dei pluviometri. Un variogramma come quello di figura, non avendo un limite superiore, indica chiaramente la non esistenza della funzione covarianza, ma, mancando l’andamento parabolico, indica anche che la media della variabile non ha un trend, senza però essere necessariamente costante nel campo. Il modello che in maniera specifica interpreta questa situazione è quello strettamente intrinseco, che, si ricorda, è un modello di FA in cui gli incrementi a distanza h sono stazionari.

Fig. 3.25 - Bacino di Kadjemeur: localizzazione dei pluviomatri e variogramma sperimentale dei mm di pioggia dopo un evento (da Mining Geostatistics di Journel & Huijsbregts, pag. 19).

Per rendersi conto di ciò riprendiamo l’andamento della FAI mostrata nella figura 3.11. Il suo variogramma (fig 3.26), calcolato per 50 passi sui sei mila dati, è lineare e non raggiunge un valore di soglia, neanche se calcolato per 5000 passi. In questa situazione, come ricordato sopra, non esiste la funzione covarianza in quanto questa, per una FAI strettamente intrinseca, non è invariante per traslazione ma dipende dalla posizione rispetto al campo. Nella figura 3.27 è riportata la covarianza (centrata) sperimentale calcolata per 50 passi sull’inera area. Ma se la stessa funzione viene calcolata in zone diverse si ottengono risultati diversi, come mostra la fig 3.28 dove sono riportate le covarianze sperimentali calcolate, per 50 passi, sui primi 2000 valori del campo, sui valori da 2001 a 4000 e sugli ultimi 2000 valori. Le covarianze sono tutte e tre diverse tra di loro e da quella calcolata su tutto il campo, riportata anche nello stesso grafico. Non così è per i variogrammi sperimentali che, calcolati sugli stessi tratti, hanno andamento molto simile tra loro e rispetto al variogramma globale (fig 3.29).

34

Variogramma

50 40 30 20 10 0 0

10

20

30

40

50

Tempo

Covarianza

Fig. 3.26 - Variogramma sperimentale delle progressive calcolato su tutti i valori della simulazione

690 685

680

675 670

665 660 0

10

20

30

40

50

Tempo Fig. 3.27 - Covarianza sperimentale delle progressive calcolata su tutti i valori della simulazione.

35

Covarianza

750

Zona 1-2000 Zona 2001-4000 Zona 4001-6000 Totale

550

350

150 0

5

10

15

20

25

30

35

40

45

50

Tempo

Variogramma

Fig. 3.28 - Covarianze sperimentali calcolate su tutto il campo e sui tratti 1-2000, 2001-4000 e 4001-6000 della simulazione

50

Zona 1-2000

40

Zona 2001-4000 Zona 4001-6000

30

Totale

20 10 0 0

10

20

30

40

50

Tempo Fig. 3.29 - Variogrammi sperimentali calcolati su tutto il campo e sui tratti 1-2000, 2001-4000 e 40016000 della simulazione

36

Un altro tipo di variogramma espressione dello stesso modello è il variogramma ad andamento logaritmico; anch’esso non è limitato superiormente ed è lontano da un andamento parabolico. 3.4.3 Comportamento nell’origine del variogramma

La continuità e la regolarità spaziale della VR sono responsabili del comportamento del variogramma nell’origine (cioè quando h → 0) . Nell’esaminare questa proprietà facciamo riferimento a variogrammi sperimentali riconducibili a FA stazionarie, quasistazionarie e strettamente intrinseche. Si possono individuare tre forme di comportamento: parabolico1, lineare, discontinuo. Nella fig. 3.30 sono riportati tre coppie di grafici: a, b, c; ogni coppia rappresenta nel grafico di sinistra l’andamento di una variabile v(x) lungo una progressiva x e nel grafico di destra il relativo variogramma.

Fig. 3.30 -Comportamento del variogramma nell’origine: discontinuo (a); lineare (b); parabolico (c) (da “An introduction to Applied Geostatistics” di E.H. Isaaks & R.M.Srivasatava). 1

Da non confondere il comportamento parabolico (o lineare) nell’origine del variogramma con il comportamento parabolico (o lineare) del variogramma nel suo complesso. Un variogramma può avere un comportamento lineare o parabolico nell’origine e nello stesso tempo ammettere un valore di soglia finito .

37

Tutti e tre i variogrammi raggiungono un valore di soglia che si mantiene poi costante, ma il comportamento all’origine è diverso. In particolare esaminando i grafici dal basso verso l’alto: • il variogramma (c) ha un comportamento parabolico nell’origine che interpreta l’elevata regolarità e continuità spaziale della variabile v di sinistra. La funzione variogramma è due volte derivabile nell’origine ed anche v(x) è derivabile; • il variogramma (b) ha nell’origine un comportamento lineare, ma non è più derivabile. La variabile v(x) è continua, ma anch’essa non derivabile; • il variogramma (a) presenta una discontinuità nell’origine: γ ( h ) non tende a zero quando h tende a zero, pur dovendo sempre essere per definizione γ ( 0) = 0 . Anche l’andamento di v(x) non è continuo, come mostra il corrispondente grafico di sinistra: la variabilità tra due valori v(x) e v(x+h) presi in due punti molto vicini può essere piuttosto elevata e questo tanto più quanto più è elevata la discontinuità nell’origine. Questa variabilità locale può essere paragonata al concetto fisico di rumore di fondo. Man mano che la distanza h aumenta, la variabilità spesso diventa più continua e ciò si riflette nella continuità di γ ( h ) per h maggiore di zero. Questa discontinuità del variogramma nell’origine viene chiamata effetto pepita (nugget effect nella terminologia inglese e effet de pépite nella terminologia francese). Nel punto successivo verrà data una spiegazione a questo comportamento.

Sempre con riferimento alla fig. 3.30 si osservi che i tre grafici v(x), campionati in maniera regolare ad un certo intervallo (sette punti campionati lungo il tratto dell’asse x mostrato), danno luogo a valori uguali (punti più marcati del grafico), pur mostrando i tre andamenti una variabilità sostanzialmente diversa. Questo fatto mette in luce un fatto molto importante: per evidenziare il comportamento all’origine del variogramma è necessario campionare il fenomeno a piccola scala. Ciò nella pratica si realizza intensificando localmente la campionatura in alcune zone, per es. secondo lo schema di fig. 3.31, dove su una campionatura a maglia regolare quadrata sono inpostate due croci di campionatura più fitta.

Fig. 3.31 - Schema di una campionatura per il riconoscimento della piccola scala.

38

La variabilità esaminata nei grafici della fig. 3.30 si riferisce ad una ricostruzione numerica (simulazione) realizzata per una rappresentazione didattica del comportamento all’origine del variogramma. Esaminiamo ora lo stesso problema su variabili realmente misurate: si tratta di uno studio realizzato da J. P. Delhomme (1976) su dati piezometrici della zona di Korhogo (Costa d’Avorio) riportato su Mining Geostatics di A. Jornail e Ch. Huijbegts. I dati si riferiscono alla profondità della falda misurata giornalmente su un certo numero di piezometri distanti tra di loro circa 500 m (vedi fig. 3.32). Per ogni piezometro era disponibile il profilo delle altezze piezometriche giornaliere per la durata della stagione delle piogge (da Agosto a Novembre).

Fig. 3.32 - Schema dell’indagine sulle piezometriche a Korhogo (da Mining Geostatics di A. Journel e Ch. Huijbegts).

La fig. 3.33 mostra quattro di questi profili relativi ai piezometri N° 3, 4, 33, 18, localizzati a profondità crescente della falda. La stessa figura nella parte bassa mostra i quattro variogrammi corrispondenti. Premesso che la quantità d’acqua trattenuta dal sottosuolo è proporzionale allo spessore di terreno attraverso cui l’acqua filtra prima di raggiungere la falda, sui profili si possono fare le seguenti osservazioni. Piezometro N° 3: la falda è molto prossima alla superficie ed il livello reagisce ad ogni precipitazione. Il profilo è molto erratico ed è direttamente legato alla precipitazione giornaliera. Il variogramma è lineare con la presenza di un effetto pepita. Piezometro N° 4: La falda è un po’ più profonda e l’effetto dell’assorbimento del sottosuolo comincia a farsi sentire. Il variogramma mostra un effetto pepita più piccolo e si intravede un valore di soglia. Piezometro N° 33: La falda è ancora più profonda e l’effetto di regionalizzazione del terreno è totale. Rimangono ancora alcune discontinuità del profilo dovute a bruschi riempimenti della falda. Il variogramma presenta ancora un piccolo effetto pepita, che

39

rappresenta i riempimenti discontinui e seguita con un andamento parabolico caratteristico di un’elevata continuità della variabile. Piezometro N. 18: La falda è ancora più profonda e gli effetti dei riempimenti discontinui sono scomparsi. E’ quindi scomparso l’effetto pepita ed il variogramma presenta un comportamento parabolico perfetto. Si fa notare che nella fig. 3.33 è mostrato solo il comportamento all’origine del variogramma e non il suo andamento globale. Ciò si deduce dal fatto che i variogrammi sono tracciati per distanze fino a 15 giorni, quando il profilo si sviluppa per 180 giorni.

Fig. 3.33 - Profili dell’altezza piezometrica e variogrammi corrispondenti (da Mining Geostatics di A. Journel e Ch. Huijbegts).

40

Un altro esempio di corrispondenza tra andamento del profilo e comportamento del variogramma è mostrato nella fig. 3.34. Essa si riferisce a due profili di tensione acquasuolo misurati uno ad una profondità di 30 cm e due giorni dell’evento irriguo, l’altro a profondità di 120 cm e dopo 14 giorni dell’evento irriguo. I due profili differiscono per la presenza in quello (B) di alcune marcate discontinuità, che sono responsabili del forte effetto pepita che compare nel variogramma (B). Il resto del variogramma ha un andamento analogo, come analoghi sono i profili, non considerando le discontinuità.

Fig. 3.34 - Profili di tensione acqua-suolo e variogrammi corrispondenti (da “spatiall variability of soilwater properties in irrigated soils”, di G. Uehara et al., Proceedings of the workshop on soil spatial variability of the ISSS and the SSSA, Las Vegas 1984).

3.4.4 Andamento e modelli del variogramma elementare

Il comportamento nell’origine di un variogramma, così come è stato descritto nel punto precedente, può essere: parabolico, lineare, discontinuo e con l’aumentare di h il variogramma aumenta di valore ed evolve sostanzialmente secondo due forme: • raggiunge un valore di soglia; • aumenta indefinitamente. Nel primo caso la FA rappresentata dal variogramma è stazionaria, o quasi- stazionaria, mentre nel secondo caso la FA è intrinseca. Il valore di soglia e la distanza alla quale esso è raggiunto vengono chiamati nella letteratura geostatistica rispettivamente: paliér e portée secondo la terminologia

41

francese e sill e range secondo la terminologia inglese (v. fig. 3.35); il relativo variogramma è perciò detto con paliér o con sill. Senza paliér o senza sill è detto invece il variogramma di una FA intrinseca. Nel seguito di queste dispense noi faremo riferimento alla terminologia inglese.

g (h) Sill

Fig. 3.35 - Variogramma con valore di soglia e significato dei parametri sill/paliér e range/portée.

Range

h

Sia che si tratti di un variogramma con sill che senza sill il suo andamento in funzione della distanza può seguire diversi stili in dipendenza del tipo di variabilità spaziale del fenomeno di studio. Piuttosto che descrivere direttamente il comportamento dei variogrammi sperimentali, cosa che sarebbe oltre che difficile poco praticabile, lo faremo attraverso l’esame delle funzioni analitiche , che, comunemente vengono assunte per descrivere il comportamento dei variogrammi sperimentali. E’ però prima necessario richiamare ed integrare le proprietà matematiche della funzione variogramma: • è una funzione che assume sempre valori positivi:

γ (h) ≥ 0 • per h=0 si ha sempre:

γ ( 0) = 0 • il variogramma è una funzione pari:

γ (h) = γ ( − h) • quando la FA è stazionaria e quindi esiste la funzione di covarianza C(h):

γ ( h ) = C ( 0) − C ( h )

(3.16)

• è dimostrato che il variogramma cresce all’infinito meno rapidamente che h 2 , cioè: γ (h) lim 2 = 0 per h → 0 h

42

• la funzione variogramma deve essere tale da dar luogo a combinazioni lineari autorizzate. Rendiamo esplicita questa condizione: Le combinazioni lineari giocano un ruolo molto importante in Geostatistica: esse possono rappresentare uno stimatore, una varianza o altre entità, e pertanto è necessario che esse siano autorizzate, cioè ammettano una varianza finita. Sia Z(x) una FA stazionaria, di media m, covarianza C(h) e variogramma γ ( h ) e sia Y una combinazione lineare finita del tipo:

Y=

n

λiZ ( xi )

(3.17)

i =1

con λi reali qualsiasi. La varianza di Y non può essere negativa: Var (Y ) ≥ 0 . Esplicitamente essa è data da: Var{Y } =

i j

λiλjC ( xi − xj ) ≥ 0 .

(3.18)

Infatti, essendo: E{Y } = E

{ }

λiZi =

i

{ 2}= i

EY

j

i

λiE{Zi} = m λi ; i

λiλjE{Z ( xi ) ⋅ Z ( xj )} ;

{ 2 }− E{Y }⋅ E{Y },

Var{Y } = E Y

e tenendo conto di: 2 C ( xi − xj ) = E{Z ( xi ) ⋅ Z ( xj )} − m , si ottiene la (3.18). La funzione C(h) deve essere tale da assicurare che l’espressione (3.18) sia non negativa. Quindi, per definizione C(h), deve essere semi-definita positiva. Tenendo conto della relazione (3.16) che lega covarianza e variogramma, la precedente può scriversi in funzione del variogramma: Var{Y } = C (0) λi λj − i j

λiλjγ ( xi − xj ) i

(3.19)

j

Se Z(x) non ammette covarianza ma solo variogramma, cioè nel caso che la FA sia intrinseca, la varianza di Y, in base alla espressione precedente esiste solo se i

λi = 0

(3.20)

ed è data da: Var{Y } = −

λiλjγ ( xi − xj ) i

j

43

Essa assume valori positivi o nulli solo se la funzione - γ ( h ) è semidefinita positiva con la condizione (3.20). Con riferimento a quest’ultima si dice anche che - γ ( h ) deve essere una funzione semidefinita positiva condizionale.

Tra le funzioni matematiche che hanno le proprietà sopra enunciate la letteratura geostatistica ne propone alcune che sin qui si sono rivelate adatte a descrivere l’andamento dei variogrammi sperimentali. Diamo di seguito le loro espressioni analitiche, classificandoli tra quelle con sill e quelle senza sill. Queste funzioni di variogramma sono anche chiamate modelli di variogramma o schemi di variogramma.

3.4.4.1 Variogrammi con sill

1. Modello pepitico

γ (h) = c 1 − δ (r )

r= h ≥0

δ(r) è una funzione che vale 1. quando r = 0 e 0. per ogni altro valore di r. Il parametro c è il sill del variogramma, che è caratterizzato dall’avere un range nullo. Questo modello esprime una discontinuità nell’origine.

2. Modello sferico

γ ( h) =

[

c 1.5 ⋅ r a − 0.5 ⋅ r 3 a 3 c

]

r= h ≤a r= h ≥a

a e c sono i parametri del modello e rappresentano rispettivamente il range ed il sill. Il comportamento nell’origine è lineare con una pendenza pari a 1.5 c/a. 3. Modello esponenziale

γ ( h) = c 1. − exp( − r a )

r= h ≥0

Il parametro c rappresenta il sill, che è per questo modello è raggiunto asintoticamente. Il modello è pertanto a range infinito. Comunque, per avere una misura della distanza entro cui si manifesta la correlazione ed anche per confronto con altri modelli è stato introdotto un range pratico a’, definito come la distanza alla quale viene raggiunto il 95% del sill. Esso risulta essere: a’= 3a. Il comportamento nell’origine è lineare con una pendenza pari a c/a.

44

4. Modello gaussiano

γ ( h) = c 1. − exp( − r 2 a 2 )

r= h ≥0

Il parametro c rappresenta il sill, che anche per questo modello è raggiunto asintoticamente con range infinito. Anche per il modello gaussiano è stato introdotto un range pratico a’ avente lo stesso significato di quello esponenziale, cioè la distanza alla quale viene raggiunto il 95% del sill. Esso risulta essere: a '= a 3 . Il comportamento nell’origine è parabolico, quindi con tangente orizzontale.

5. Modello a effetto buco

γ (h) = c[1 − exp(− r a) ⋅ cos(2πfr )] r = h ≥ 0 Questo modello ha l’andamento di un esponenziale, ma con una forma oscillatoria di periodo 1/f. I parametri a e c hanno lo stesso significato del modello esponenziale. Affinchè questa funzione abbia le caratteristiche matematiche sopra enunciate è necessario che in R2 f e a soddisfino il seguente vincolo: 1 f ≥ 2 π a . Questa restrizione non è richiesta in R1.

3.4.4.2 Variogrammi senza sill

Questi modelli, si ricorda, corrispondono a FA che hanno una illimitata capacità di dispersione e quindi non ammettono funzione covarianza. I modelli più noti sono quelli appartenenti alla classe seguente: 1. Modelli potenza

γ (h) = r θ ; r = h > 0 e 0 < θ < 2 Quello più correntemente usato nella pratica è il modello lineare:

γ ( h ) = ωr r = h > 0 , dove ω è la pendenza nell’origine. Man mano che nel modello potenza θ aumenta da 1 a 2, il comportamento di γ ( h ) diventa sempre più regolare nell’origine, rappresentando una FA sempre più regolare. In questa situazione però spesso non si riesce a distinguere se si tratta di una FA stazionaria che ammette un variogramma potenza con valore di θ prossimo a 2, o se si tratta di una FA con deriva. La scelta è lasciata alla sensibilità dell’utente.

45

3.4.4.3 Altri tipi di variogramma

Vi è un modello di variogramma che gioca un importante ruolo nello studio geostatistico dei fenomeni periodici. Questo modello si chiama appunto periodico, non rientra nelle due categorie precedenti e viene presentato a parte. 1. Modello Periodico

γ ( h ) = ξ 1 − cos( 2 π T ⋅ r ) ; ξ ≥ 0; r = h ≥ 0; E’ un variogramma che non si stabilizza su un valore di soglia, nemmeno all’infinito, quindi non ammette sill e la correlazione spaziale non si estingue mai. Nello stesso tempo non cresce indefinitamente, ma è limitato ed il limite superiore vale ξ. Tuttavia è possibile associare al variogramma una varianza pari a ξ. Il parametro T è il periodo del variogramma. Per ognuno dei modelli di variogramma presentati, nei grafici seguenti viene mostrato l’andamento della funzione analitica e l’andamento di una variabile spaziale che può essere descritta da una FA avente come funzione strutturale quel variogramma. Nelle figure, le funzioni variogramma sono riportate a destra e le variabili spaziali a sinistra.

46

Fig. 3.36 – Variogrammi modello

47

Fig. 3.37 – Variogrammi modello

48

3.4.5 Variogrammi a strutture annidate

E’ frequente nella pratica delle applicazioni riscontrare nell’andamento di un variogramma sperimentale variazioni di pendenza. Per rendersi conto si possono osservare i variogrammi di fig. 3.38. Essi rappresentano, dall’alto verso il basso, i variogrammi del contenuto in Ca, Cl e NO3 disciolti in acqua misurati in 68 sorgenti del bacino di Dyle (Belgio), in tre anni diversi:1975, 1981, 1983. Senza ora entrare nel merito dell’interpretazione (lo faremo un po’ più avanti) ci limitiamo ad osservare che tutti i variogrammi presentano un cambiamento di pendenza, per alcuni più accentuato per altri meno, che si verifica sistematicamente attorno al chilometro. Un altro cambiamento di pendenza ha luogo, per alcuni variogrammi più visibile per altri meno visibile, attorno agli otto chilometri, quando viene raggiunto il valore di soglia. Inoltre alcuni dei variogrammi mostrano un debolissimo effetto pepita. Questo comportamento, che è riconducibile alla particolare fenomenologia, ha una spiegazione molto semplice.

Fig. 3.38 - Variogrammi sperimentali dei contenuti in Ca, Cl e NO3 misurati in tre diversi periodi (da “study of spatial and temporal variations of hydrogeochemical variables using factorial kriging analysis” di P. Goovaerts & Ph. Sonnet, Proc. 4th Int. Geostat. Congress, Troia, 1992).

49

Il cambiamento di pendenza è un segno che il variogramma sperimentale è costituito dalla sovrapposizione di più variogrammi elementari, aventi diversi valori di range. In particolare i cambiamenti di pendenza si verificano in corrispondenza dei ranges dei variogrammi elementari. Ciò appare evidente osservando gli schemi di fig. 3.39. Là dove i variogrammi di fig. 3.38 non presentano effetto pepita essi possono quindi essere considerati come la somma di due variogrammi, uno con range 1 Km e l’altro con range 8 Km. Per i variogrammi che presentano effetto pepita ai due variogrammi elementari precedenti bisogna aggiungere un altro variogramma elementare: un effetto di pepita puro. E’ evidente che i variogrammi elementari perchè possano essere riconosciuti dai cambiamenti di pendenza devono essere dei variogrammi che ammettono sill, perchè solo essi ammettono range. I variogrammi elementari che compongono variogrammi di questo tipo vengono anche chiamati strutture spaziali, ognuna rappresentante di una scala spaziale paragonabile al proprio range. Poiché le scale spaziali sono progressive, queste strutture sono chiamate annidate (nested nella terminologia inglese e gigognes in quella francese) ad indicare che esse sono inscatolate l’una dentro l’altra. Anche i variogrammi risultanti da strutture annidate vengono chiamati annidati. E’ chiaro che nel corso di una indagine le strutture che possono essere evidenziate sono quelle che si esprimono alle scale del lavoro, che vanno da distanze pari al lato della maglia di campionatura a distanze pari alle dimensioni e forse anche più del dominio di lavoro: le strutture a distanze più elevate non si colgono, mentre le strutture a distanze più piccole si confondono.

g 1(h)+g 2(h)

g 1(h)

g 2(h)

h

h

g 1(h)+g 2(h)+g 3(h) g 1(h) g 2(h) g 3(h) h

h

Fig. 3.39- Sovrapposizione di variogrammi elementari.

Ora, poichè un variogramma elementare è l’espressione di un fenomeno particolare a cui è da attribuirsi l’esistenza e la strutturazione spaziale di una variabile, un variogramma annidato può essere considerato come l’espressione di più fenomeni

50

sovrapposti, anzi, dimostreremo quì di seguito, che una variabile che presenta un variogramma annidato può essere considerata come la somma di più variabili indipendenti. Sia Z(x) una FA stazionaria che ammette variogramma γ ( h ) . Supponiamo che essa sia composta dalla somma di S variabili indipendenti, ognuna Z u ( x ) stazionaria di variogramma γu ( h ) . Si ha: Z ( x) =

S −1 u=0

(3.21)

Zu( x)

e anche: Z ( x + h) =

S −1 u =0

Zu ( x + h) .

Il variogramma di Z(x) è data pertanto dalla seguente espressione:

{

}

γ (h) = E [Z ( x + h) − Z ( x)]2 = E

S −1 u =0

(Zu ( x + h) − Zu( x))

2

=

S − 1 S −1 u = 0 u '= 0

E{[Zu ( x + h) − Zu ( x)] ⋅ [Zu '( x + h) − Zu '( x)]}

Separando nella precedente i termini rettangolari si ha che:

{

}

E [Zu ( x + h) − Zu ( x)] +

γ (h) =

2

u≠

E{[Zu ( x + h) − Zu ( x)] ⋅ [Zu '( x + h) − Zu '( x)]} u'

Ora essendo le variabili Zu(x) spazialmente indipendenti i termini rettangoli sono nulli e quindi si ottiene:

γ (h) =

S −1 u =0

{

}

E [Zu ( x + h) − Zu ( x)] = 2

S −1

γu (h) ,

u =0

cioè il variogramma di Z(x) è costituito dalla somma dei variogrammi delle variabili Zu(x). Queste ultime variabili vengono anche chiamate componenti spaziali di Z(x), perchè ognuna si esprime ad una determinata scala spaziale. Il sill di ogni componente spaziale rappresenta la quota parte della dispersione della variabile che compete a quella scala. Alla luce di questo fatto possiamo dare un’interpretazione di ciò che è stato chiamato effetto pepita. Supponiamo che una variabile, il cui variogramma mostra un effetto pepita, abbia un certo numero di strutture spaziali a scale inferiori del lato della maglia, o inferiori alla più piccola distanza tra i campioni, e che ognuna di queste strutture spaziali sia rappresentata da un variogramma con sill. Sommando tutti questi variogrammi e rappresentando la somma alla scala dei variogrammi sperimentali, si ha che essa è visibile come una discontinuità nell’origine di entità pari alla somma dei sills. Quindi quando si riscontra in un variogramma sperimentale con una discontinuità nell’origine, questa può essere dovuta alla variabilità espressa da un numero imprecisato

51

di strutture spaziali presenti a scale più piccoli di quelle di lavoro. Queste strutture sono confuse in un unico rumore di fondo. Un’altra causa che concorre alla formazione di un effetto pepita sono gli errori di misura. L’errore di misura può infatti, anche se non ha carattere spaziale, essere considerato come una componente della variabile bruta: una delle componenti Z u '( x ) di Z(x) secondo l’espressione (3.21). Poichè Z u '( x ) è indipendente da ogni altra variabile, il variogramma di Z(x) comprende il variogramma degli errori di misura, che, essendo indipendenti spazialmente, presentano come variogramma un effetto di pepita puro a tutte le scale. Siamo ora in grado di dare una interpretazione fenomenologica alla variabilità espressa dai variogrammi sperimentali di fig. 3.38. Tale interpretazione è quella contenuta nel lavoro citato. La struttura a piccola scala può essere dovuta alla presenza di sorgenti locali di inquinamento come aziende agricole, discariche di rifiuti solidi o aree urbanizzate, mentre la struttura a grande scala è principalmente da attribuirsi alla variabilità delle caratteristiche geologiche dell’acquifero che ha luogo su vasta area. I contributi delle due strutture ai variogrammi riflettono il ruolo giocato dai fattori geologici e antropici nel comportamento spaziale di ogni ione. Infatti mentre la provenienza di NO3 e Cl può essere attribuita ad attività umane, il contenuto di Ca nelle sorgenti d’acqua è controllato su vasta scala dalla presenza nell’acquifero di carbonati. Il comportamento dei variogrammi mostra inoltre che la componente a grande scala del Ca si mantiene inalterata nel tempo, mentre quella a piccola scala si indebolisce. Quanto al Cl e NO3 il contributo della granda scala varia sensibilmente con gli anni, mentre quello a piccola scala si mantiene sostanzialmente invariato. Diamo ancora altri esempi di variogrammi annidati: quello di fig. 3.40 si riferisce al variogramma del Co, costruito sui dati trasformati (trasformazione gaussiana) provenienti da una campionatura dei suoli del sud-est della Scozia. La campionatura, che ha interessato anche il Cu, è stata condotta per lo studio della deficienza dei due metalli della regione scozzese. Come si può osservare il variogramma presenta un forte effetto pepita, una struttura locale a 3 km e mezzo ed un’altra a 16 km. Della dispersione totale il 60% si consuma a piccola scala, un altro 25% alla scala dei 3 km e mezzo ed il resto alla scala dei 16 km.

52

Fig. 3.40 - Variogramma del Co nei suoli del sud-est della Scozia (da “Disjunctive kriging in agricolture” di R. Webster and M.A. Oliver in M. Armstrong (ed), Geostatistisc, vol.I, pag.427).

I variogrammi di fig. 3.41 si riferiscono ai contenuti in Cu, Pb e Zn di campioni di suolo provenienti da una campagna di prospezione geochimica nella regione di Munster (Francia). I variogrammi mostrano la presenza di tre strutture spaziali per tutte e tre le variabili. Per il Pb: una pepitica, una a 1.5 km e l’altra a 10 km. Per lo Zn: una pepitica, una a 2.0 km e l’altra a 7.5 km. Per il Cu: una pepitica, una a 2.0 km e l’altra a 9.0.

53

Fig. 3.41 - Variogrammi di Cu, Pb e Zn di una campagna geochimica nella regione di Munster (Francia) (da “Analyse krigeante de données geoquimiques” di L. Sandjivy in Sciences de la Terre, n.24).

3.4.6 Anisotropie nei variogrammi Le VR sono quasi sempre riferibili a fenomeni naturali od a fenomeni in cui è presente una componente naturale. Tali fenomeni spesso si sviluppano in maniera differente in funzione della direzione: un fenomeno quasi sempre presenta delle direzioni principali. Lo si è già visto nel paragr. 3.4.2 dove in alcuni esempi mostrati il fenomeno si presentava non stazionario in una direzione e stazionario in quella ortogonale. Si tratta di una anisotropia che riguarda la presenza di deriva, che quasi sempre per sua natura si svolge secondo una determinata direzione. Le anisotropie però di cui ci occupiamo in questo paragrafo riguardano il comportamento del variogramma in situazioni di stazionarietà. Esso infatti può avere lo stesso andamento in tutte le direzioni, ed allora la funzione variogramma è isotropa. Si osservi per esempio la fig. 3.42 in cui sono mostrati i variogrammi sperimentali della riflettanza di una scena Landsat TM7 (Path 182,row 066 del 14.06.2000). I variogrammi si riferiscono alle quattro direzioni geografiche principali. In questo caso il variogramma medio, vale a dire quello calcolato lungo una (qualsiasi) direzione con una tolleranza angolare di 90°, rappresenta la funzione variogramma in tutto il piano.

Fig. 3.42 – Variogrammi secondo le quattro direzioni principali (N-S, NE-SW, NE-SW) della riflettanza di una porzione di scena Landsat TM7 (Path 182,row 066 del 14.06.2000).

Si osservi invece la fig. 3.43. Essa mostra i variogrammi sperimentali relativi alle variabili Cr, Cu e Ni misurate su campioni di sedimenti di lago dell’area di Schefferville

54

(Quebec Nord). Come si può notare il sill è uguale a 3.2 per entrambe le direzioni, mentre il range è 7 km per la direzione N58E (B) e 13 km per la direzione N32W (A). Questa differenza di comportamento che riguarda il range lasciando il sill invariato è noto in geostatistica come anisotropia geometrica.

Fig. 3.43 - Anisotropia geometrica del variogramma multivariabile calcolato sui dati di Schefferville (fig. tratta da “Spatial filtering under the linear coregionalization model” di G. Bourgault e D. Marcotte, Troia 92).

Un altro esempio di anisotropia geometrica si può osservare nella fig. 3.44, che riporta i variogrammi orizzontale e verticale della conducibilità idraulica derivati dall’analisi granulometrica di campioni raccolti in un sito della Germania del Nord ad una distanza orizzontale di 2.5 m e verticale di 0.5 m. I valori dei ranges che si osservano sono 1.5 e 14 m. Si osservi ora la fig. 3.45. Essa si riferisce ai variogrammi relativi al contenuto in metallo Pb+Zn di una miniera italiana nelle Alpi Orientali. Lo schema di sinistra riporta le direzioni indagate. Le direzioni 3 e 4 avendo dato dei variogrammi con andamento simile sono stati mediate ed il variogramma medio ottenuto assieme a quello della direzione 1 sono mostrati nella figura. Vi si nota una marcata differenza di comportamento tra i due variogrammi mostrati. Quello relativo alla direzione 1 presenta un sill molto più elevato

Fig. 3.44 - Variogrammi orizzontale e verticale della conducibilità idraulica manifestanti una anisotropia geometrica (da “Spatial structure of hidraulic conductivity in various porous mediaproblems and experien-ces” di M.Th. Schafmeister & A. Pekdeger, Troia, 92)

55

dell’altro variogramma, che corrisponde ad una direzione suborizzontale. Una tale differenza nei variogrammi, che riguarda il grado di variabilità con la direzione, prende in geostatistica il nome di anisotropia zonale. Essa è caratterizzata dall’avere una direzione in cui la variabilità è massima e tale direzione è chiamata appunto direzione di zonalità. Nel caso di figura la direzione di zonalità è la direzione 1, che corrisponde alla direzione ortogonale alla giacitura di stratificazione dei calcari in cui la mineralizzazione è insediata.

Fig. 3.45 - Esempio di anisotropia zonale in una miniera di Pb e Zn delle Alpi Orientali (fig. tratta da “evaluation and optimization of a metal mine” di M. Guarascio e G. Raspa, Proceedings Apcom 74, Denver).

Un altro esempio di anisotropia zonale è riportato nella fig. 3.46. Esso si riferisce ad uno studio dei suoli in un bacino della Costa d’Avorio. La variabile studiata è l’impoverimento di sostanze umiche nel terreno, misurata dal rapporto tra la quantità di argille e limi presenti negli orizzonti d’accumulazione e quella presente negli orizzonti umici. La parte superiore della figura mostra lo schema di campionatura; esso è costituito da una serie di allineamenti subparalleli ortogonali all’asse del bacino, che è orientato NE-SW. I grafici B e C sottostanti si riferiscono ai variogrammi lungo la direzione NE-SW e lungo la direzione ortogonale. Il variogramma NE-SW mostra una struttura elementare caratterizzata da un range di 300-400 m ed un sill di 0.5, mentre Il variogramma della direzione NW-SE mostra chiaramente la sovrapposizione di due strutture una di range attorno ai 300 e sill 0.5 (la stessa che compare nel variogramma precedente), l’altra di range 1240 e sill 0.35.Si può pertanto affermare che il fenomeno presenta due strutture spaziali: una isotropa di range 300 m e l’altra anisotropa zonale di range 1240 m e direzione di zonalità NW-SE.

56

Quest’ultima direzione, essendo ortogonale all’asse del bacino, si conferma essere appunto quella di maggiore variabilità. Nel paragrafo 3.4.6 sarà mostrato come procedere alla modellizzare delle anisotropie della funzione variogramma.

Fig. 3.46 - Schema di campionatura e variogrammi della carenza in sostanze umiche in un bacino della Costa d’Avorio (da “Analyse de la variabilité d’un parametre pedologique” di J.M.Iris, in Sciences de la Terre, Vol 24, 1984).

3.5 Aggiustamento di un variogramma modello Affinché quanto espresso dai variogrammi sperimentali possa essere usato, quale modello di variabilità spaziale, per svolgere le operazioni geostatistiche che sono richieste nell’ambito di un progetto, è necessario che sia trasformato in una funzione analitica γ ( h ) capace di fornire un valore del variogramma in funzione della distanza e dell’orientazione di una qualsiasi coppia di punti dello spazio. In pratica se x1 e x2 sono le posizioni di due punti in R2 (definite dalle coordinate x1u, x1v; x2u, x2v), la funzione variogramma deve potersi esprimere come funzione delle componenti: hu = x2u-x1u e hv = x2v-x1v. Inoltre la funzione - γ ( h ) deve essere semidefinita positiva, cioè tale che, per qualunque n finito e maggiore di zero e per qualunque set di numeri reali { λi} (i = 1, n) associati ad altrettanti punti dello spazio { xi} , comunque posizionati, si abbia:

57

λiλjγ ( x i − x j ) ≥ 0 ,

i

j

λi = 0 se il modello è solo intrinseco.

con la condizione i

Sulla base di quanto è stato detto nel paragrafo 3.4.5 a proposito delle strutture annidate, la funzione γ ( h ) in generale può considerarsi come costituita dalla somma di s variogrammi elementari, ognuno espresso da una delle funzioni presentate nel paragrafo 3.4.4. Utilizzeremo d’ora in poi la notazione seguente:

γ (h) =

s −1

γu ( h )

u=0

dove u è l’indice della struttura elementare e γu il suo variogramma. Nella precedente si è inteso indicare con l’indice 0 la struttura pepitica, quasi sempre presente. Se non vi è effetto pepita la struttura può essere formalmente mantenuta attribuendo ad essa un sill nullo. In queste condizioni aggiustare un modello vuol dire: ♦ definire il numero s di strutture; ♦ per ogni struttura specificare: se essa è isotropa o anisotropa se è anisotropa quale è il tipo di anisotropia (geometrica, zonale) i parametri dell’anisotropia il tipo di funzione modello (sferica, exp, gauss, etc) i valori dei parametri della funzione (sill, range) L’aggiustamento viene effettuato operando un fitting tra la funzione analitica γ ( h ) così composta ed il variogramma sperimentale. Le funzioni variogramma elementare ed il significato dei loro parametri sono stati già descritti nel paragrafo 3.4.4. Introduciamo ora i parametri che caratterizzano le anisotropie e diamo l’algoritmo per il calcolo di γ ( h ) in funzione di detti parametri. Anisotropia geometrica Questo tipo di anisotropia ricorre per di più nei fenomeni bidimensionali. In esso il range del variogramma varia con la direzione, presentandone una dove il range è massimo e l’altra, perpendicolare alla prima, dove il range è minimo. In funzione della direzione il range varia come il raggio vettore di una ellisse, i cui assi maggiore e minore, sono appunto i ranges massimo e minimo. I parametri dell’anisotropia, sono facilmente identificabili: • angolo ϕg che la direzione di massimo range forma con l’asse u del riferimento di lavoro; • rapporto λ tra range maggiore e range minore; Dato un vettore h, di componenti hu e hv, la maniera più usuale di calcolare il valore di γ(h) in funzione dei parametri precedenti è di operare come segue: 58

a) effettuare una rotazione di ampiezza ϕg degli assi coordinati fino a far coincidere l’asse u con l’asse maggiore dell’ellisse. La matrice di trasformazione è in questo caso:

R ϕg

cos( ϕ

=

− sin( ϕ

g g

)

sin( ϕ

)

cos( ϕ

g

)

g

)

Pertanto: hu'''

= [ Rϕg ] ×

' ' ' v

h

hu ; hv

b) trasformare l’ellisse in un cerchio avente raggio uguale all’asse maggiore dell’ellisse. Ciò si ottiene moltiplicando h ''' per il rapporto di anisotropia. La matrice di trasformazione corrispondente sarà: 1 0 . 0 λ

[λ] =

Pertanto: hu'' hv''

= [λ] ×

hu''' hv'''

;

c) riportare gli assi coordinati nella posizione originaria effettuando una rotazione di ampiezza -ϕg , opposta alla precedente. La matrice di trasformazione sarà: R − ϕg

=

cos(ϕg ) − sin ( ϕg ) sin ( ϕg ) cos( ϕg )

e pertanto: hu' hv'

= [ R − ϕg ] ×

hu'' hv''

.

La trasformazione complessiva sarà: hu' hv'

= [ R − ϕg ] × [ λ ] × [ Rϕg ] ×

hu hu = [ A] × hv hv

59

con A=

a c c b

e a = cos2 ϕg + λsin2ϕg b = sin2ϕg + λ cos2 ϕg c = (1 − λ ) sin ( ϕg ) ⋅ cos ϕg

Ottenuti hu' e hv' , il valore di γ(h) = γ(hu,hv) si calcola utilizzando una funzione isotropa con range pari a quello massimo e con valore di h trasformato: 2

2

γ (hu , hv ) = γ ' hu' + hv' .

Anisotropia zonale Qualsiasi anisotropia non riconducibile ad una trasformazione lineare delle coordinate può essere ricondotta a questo modello di anisotropia. Esso prevede che la variabilità si manifesti secondo una particolare direzione, che è detta direzione di zonalità. Questa direzione, individuata dal versore z, è, a due dimensioni, caratterizzata dall’angolo ϕz che essa forma con l’asse u. Nella anisotropia zonale la funzione variogramma si esprime analiticamente nel modo seguente:

γ ( h ) = γ ( hz ) , dove hz è la componente del vettore h nella direzione di zonalità: hz = h⋅z = |h|⋅⋅cos(ϕz). Si è ora in possesso di tutti gli elementi per calcolare il valore analitico di γ(h) per qualsiasi tipo di variogramma. L’aggiustamento del modello complessivo può essere fatto in due modi: Interattivamente Attraverso un opportuno programma di calcolo, vengono forniti dall’utente i parametri del modello, cioè tutte le grandezze sopra specificate. Il programma calcola i valori di γ(h) nelle direzioni indicate e l’utente, sempre attraverso lo stesso programma, li confronta con i valori sperimentali. Procedendo per tentativi l’utente modifica interattivamente i parametri del modello fino a quando non è soddisfatto del fitting, valutato visivamente (al video).

60

Automaticamente L’aggiustamento può essere automatizzato nella definizione dei sills. L’utente fornisce interattivamente il numero di strutture e, per ogni struttura: la funzione modello ed il relativo range; il tipo di anisotropia (se presente) ed i suoi parametri. I sills delle strutture possono essere calcolate minimizzando con il criterio dei minimi quadrati gli scarti tra i valori calcolati e quelli analitici del modello. La quantità da minimizzare può essere definita come segue: siano • nd il numero di direzioni (id=1,nd) su cui sono stati calcolati i variogrammi sperimentali o che si vogliono far intervenire nell’aggiustamento; • pas(id) il passo di calcolo del variogramma per la direzione id; • np(id) il numero di passi per la direzione id su cui è stato calcolato il variogramma o che si vogliono considerare ai fini dell’aggiustamento. • γ id* ( k ⋅ pas(id )) il variogramma sperimentale a passo k della direzione id; • γ id ( k ⋅ pas(id )) il variogramma analitico (calcolato con il modello) a passo k della direzione id. Esso è funzione dei sills.

La quantità da minimizzare è: nd np ( id ) id =1 k =1

[

]

W ( k , id ) γ id* ( k ⋅ pas(id ) ) − γ id ( k ⋅ pas(id ))

2

dove W(k,id) è un opportuno coefficiente per tener conto del numero delle coppie del variogramma sperimentale e dei primi passi rispetto agli ultimi. Data la complessità della espressione precedente, la procedura di minimizzazione conviene che sia iterativa, piuttosto che regressiva. Un esempio di aggiustamento lo si può avere considerando i variogrammi delle figure 3.46b e 3.46c. In questo esempio i variogrammi sperimentali interessano due sole direzioni: NE-SW e NW-SE. Sulla base di quanto espresso come commento ai variogrammi sperimentali, si può aggiustare un modello a due strutture: • la prima isotropa, di modello sferico di range 300 m e sill 0.5; • la seconda anisotropa zonale con direzione di zonalità NW-SE modellizzata con un variogramma sferico di range 1250 m e sill 0.5.

3.6 Regolarizzazione di una FA Nel paragrafo 3.1 è stato introdotto il concetto di supporto. Esso, ripetiamo, è l’entità geometrica, caratterizzata da forma e dimensioni, su cui è definita una variabile. Con riferimento a questo problema quanto è stato esposto nei paragrafi precedenti è del tutto generale e vale per qualsiasi supporto anche se, per semplicità di esposizione, si è fatto

61

sempre riferimento a variabili puntuali, cioè a supporto puntuale. E’ stato anche detto nel paragrafo citato e fatto osservare nella fig. 3.1 che se si passa da una VR puntuale ad un’altra di supporto più grande si ottiene una VR più regolare, ma in qualche forma legata alla prima. Questo passaggio da un supporto puntuale ad un altro di dimesioni più grandi è noto nella terminologia geostatistica con il nome di regolarizzazione e la variabile corrispondente è chiamata variabile regolarizzata. Questo paragrafo sarà dedicato alla ricerca del formalismo di passaggio riguardante i primi due momenti delle FA coinvolte. Il passaggio invece che interessa la legge monovariabile e la legge bivariabile viene trattato nella Geostatistica non lineare. Per tutte queste operazioni concernenti la regolarizzazione opereremo nel quadro dei modelli stazionari e quasi-stazionari. Sia Z(x) una FA puntuale stazionaria di media m, covarianza C(h), variogramma γ(h) e dominio di definizione S. In ogni punto x di S si può definire una nuova FA Zv(x), derivata dalla quella puntuale: Zv ( x ) =

1 Z (ξ )dξ v vx

dove vx è un elemento di volume centrato in x e v il suo volume. Chiamando con mv, Cv(h) e γv(h) rispettivamente la media, la covarianza ed il variogramma di Zv(x), vediamo di seguito come esse possono essere espresse in funzione delle equivalenti grandezze di Z(x). Media mv = E {Zv ( x )} =

1 E { Z (ξ )}dξ v vx

Poichè Z(x) è stazionaria, E {Z ( x )} = m , e quindi si ha: mv = m Covarianza Cv (h) = E [ Zv ( x + h) ⋅ Zv ( x )] − m2

(3.22)

Nel termine di destra Zv(x) e Zv(x+h) sono due variabili aleatorie definite su supporti di volume v centrati nei punti x e x+h (vedi fig. 3.47).

62

Fig. 3.47 - Supporti delle variabili Zv(x) e Zv(x+h).

Indicando con: Zv ( x ) =

1 Z (ξ )dξ v vx

(3.23)

e Zv ( x + h ) =

1 Z (ξ ')dξ ' v vx + h

e sviluppando la (3.22) si ottiene: Cv (h) =

1 v2

vx vx + h

E [ Z (ξ ) ⋅ Z (ξ ')]dξdξ ' −m2 =

1 v2

vx vx + h

C (ξ − ξ ')dξdξ ' .

(3.24)

L’ultimo termine a destra indica il valor medio che la funzione covarianza puntuale C(h) assume quando gli estremi del segmento ξ − ξ ' variano, in maniera indipendente, l’uno in vx e l’altro in vx + h (fig. 3.47). Essendo la funzione C(h) dipendente solo dalla distanza della coppia di punti ξ − ξ 'e non dalla particolare posizione, si ha, per la (3.24), che anche Cv(h) non dipende dal punto di appoggio, e quindi Zv(x) è anch’essa stazionaria. Per h=0, Cv(0) rappresenta la varianza della variabile regolarizzata Zv(x). Essa, in base alla (3.24), risulta essere: Cv (0) =

1 v2

vx vx

C (ξ − ξ ')dξdξ ' .=

1 v2

v v

C (ξ − ξ ')dξdξ '

(3.25)

Sostituendo nella precedente C( ξ − ξ ' ) con la sua espressione in funzione del variogramma si ottiene: Cv (0) = C (0) −

1 v2

γ (ξ − ξ ')dξdξ '

(3.26)

v v

Essa mostra che la varianza della variabile regolarizzata è più bassa di quella della variabile puntuale e la differenza è data dal termine: 63

1 v2

γ (ξ − ξ ')dξdξ ',

(3.27)

v v

che dipende dalla funzione variogramma e dal volume v, nella sua forma e dimensioni. Esso, entro dimensioni più piccole del range del variogramma, aumenta con le dimensioni di v e conseguentemente la varianza di Zv(x) diminuisce, così come deve essere. E’ molto ricorrente nella letteratura geostatistica usare per i termini integrali che compaiono nelle espressioni (3.24), (3.25), (3.26) una espressione più sintetica, tramite la quale la (3.24) si scriverebbe per esempio: Cv ( h ) = C ( v x , v x + h )

(3.28)

dove la barra sulla C indica che si tratta di un valor medio della funzione C(h) e gli argomenti indicano i supporti su cui gli estremi del segmento sono di volta in volta fatti variare in maniera indipendente. Data la stazionarietà di Zv(x) la (3.28) si scrive in forma più generale: Cv ( h ) = C ( v , v h )

Variogramma

γ v (h) =

{

1 2 E [ Zv ( x + h ) − Zv ( x ) ] 2

}

(3.29)

Sostituendo nella precedente le Zv con le espressioni (3.23) si ottiene:

γv (h) =

1 E 2

γv (h) =

1 1 2 v2

−

1 2 2 v2

vx vx + h

1 v

vx

vx

Z (ξ ) dξ −

vx

1 v

vx + h

2 Z (ξ ')dξ ' ;

E [ Z(ξ) ⋅ Z(ξ ' ) ]dξdξ ' +

1 v2

E [ Z(ξ) ⋅ Z(ξ ' )]dξdξ ' ;

64

vx + h

vx + h

E [ Z(ξ) ⋅ Z(ξ ' ) ]dξdξ ' +

γv (h) = −

1 2 2 v2

1 1 2 v2 vx vx + h

vx

vx

C(ξ − ξ ' )dξdξ '+

1 v2

vx + h

vx + h

C(ξ − ξ ' )dξdξ ' +

C(ξ − ξ ' )dξdξ ' ;

(3.30)

Data la stazionarietà di Z(x) i primi due termini della espressione precedente sono uguali e quindi, adottando la notazione sintetica, l’espressione (3.30) si può scrivere:

γv ( h ) = C ( v x , v x ) − C ( v x , v x + h )

(3.31)

) con C(0) - γ( ξ − ξ ' ) ed usando anche per il Sostituendo nella (3.30) C( ξ − ξ ' variogramma la notazione sintetica, la (3.31) diventa:

γv ( h ) = γ ( v x , v x + h ) − γ ( v x , v x ) e, data la stazionarietà di Zv(x):

γv ( h ) = γ ( v , v h ) − γ ( v , v )

(3.32)

Poichè la FA Z(x) è dotata di sill, per h= ∞ la precedente diventa:

γv ( ∞) = γ ( ∞) − γ ( v , v )

(3.33)

La precedente è analoga alla (3.26) e mostra come il sill del variogramma regolarizzato è più basso di quello del puntuale della quantità γ ( v , v ) . E’ la stessa quantità (cfr. la 3.25) di cui si abbassava la varianza per effetto della regolarizzazione. Non può essere altrimenti dato che, in condizioni di stazionarietà, la varianza di una FA, regolarizzata o no, rappresenta sempre il sill del variogramma. Nella operazione di regolarizzazione anche il range del variogramma subisce delle variazioni e queste dipendono esclusivamente dalla forma e dalle dimensioni del supporto di regolarizzazione. In particolare passando da una FA puntuale ad una regolarizzata di supporto v convesso, il range del variogramma in una certa direzione aumenta esattamente del diametro di v nella direzione considerata. La spiegazione è immediata se si pensa alla definizione di range (esso è dato dalla distanza più piccola alla quale due variabili spaziali sono indipendenti) e al fatto che due variabili regolarizzate (cioè due variabili definite su due volumi posti ad una certa distanza) sono indipendenti quando tutte le coppie di variabili associate a tutte le coppie di punti, appartenenti un punto ad un volume, l’altro punto all’altro volume, sono indipendenti. Questo concetto è espresso nella fig. 3.48.

65

Fig. 3.48 - Ranges e regolarizzazione.

3.6.1 Calcolo delle funzioni γ

Il calcolo delle espressioni intergrali che compaiono nel punto precedente viene nella pratica effettuato non analiticamente (perchè sarebbe molto complicato trattandosi, in R2, di integrali quadrupli di funzioni del tipo esaminate 3.4.4) ma numericamente. Ciò si realizza discretizzando i volumi in insiemi di punti e trasformando gli integrali in sommatorie. Supponiamo per es. di dover calcolare l’espressione:

γ (v1, v 2) =

1 v1v 2

v1 v 2

γ (ξ − ξ ')dξdξ

dove v1 v2 sono due volumi di forma quasiasi, diversi tra loro e comunque posizionati nello spazio. E’ questo un caso generale che comprende tutte le configurazioni sopra esaminate, dove i due volumi sono sempre uguali, posizionati alcune volte in maniera coincidente, altre volte ad una certa distanza tra di loro. Si discretizzano dapprima, con la stessa modalità, i due volumi in due insiemi di punti di numero rispettivamente n1 e n2 e se ne calcolano le coordinate. Si considerano poi tutte le coppie di punti, il primo nel volume v1, il secondo nel volume v2 e per ogni coppia si calcola la distanza ed il valore corrispondente del variogramma. In fine si mediano su tutte le coppie i valori del variogramma. Questo tipo di operazione è illustrata nella figura 3.49. La precisione nel calcolo dipende dalla densità di discretizzazione rispetto alla forma dei volumi. Per figure regolari quali blocchi rettangolari una discretizzazione in 16 punti dà, per i problemi più comuni, risultati accettabili.

66

γ (v1, v 2) =

Fig. 3.49 - Calcolo numerico dell’espressione

1 n1n 2

n1

n2

γ ( xi − xj )

i =1 j =1

γ.

3.6.2 Regolarizzazione di un effetto pepita . Supponiamo che il variogramma da regolarizzare si componga di S strutture elementari:

γ (h) =

ns −1

γu ( h )

u=0

e tra queste una microstruttura (componente pepitica) γ 0( h ) . Il valore di γv ( h ) si ottiene normalmente dalla (3.32) sommando il contributo di tutte le strutture, calcolato come indicato nel paragrafo 3.6.1. E’ però possibile ricavare il contributo della componente pepitica γov ( h ) direttamente dalla 3.32 in maniera molto immediata. Riportiamo di seguito il metodo a solo scopo didattico, essendo dal punto di vista pratico più conveniente eseguire il calcolo secondo la procedura già citata, se non altro per conformità con le altre strutture. Consideriamo nella 3.32 una distanza h più grande delle dimensioni di v, cioè tale che i volumi v e vh siano disgiunti. Applichiamo la 3.32 alla sola componente pepitica:

γ 0 v ( h ) = γ 0 ( v , v h ) − γ 0( v , v ) .

(3.34)

Per comodità operiamo con la covariaza, anzicchè con il variogramma:

γ 0 (h ) = C 0 (0) − C 0 (h ).

(3.35)

C 0 (h ) è chiaramente una funzione di covarianza con un range a0 molto piccolo rispetto alla scala del lavoro, quindi anche rispetto alla taglia di v (v.fig. 3.50). 67

C(h) C(0)

a0 << v Fig. 3.50 - Funzione Covarianza della microstruttura.

h

a

I due termini di destra della 3.34 nella loro forma esplicita e con la sostituzione data dalla 3.35 precedente, diventano:

γ 0 ( v , vh ) =

1 v2

γ 0( v , v ) =

1 v2

v v+h

v v

C 0(0) dξdξ '−

C 0(0) dξdξ '−

1 v2

1 v2

v v

v v+h

C 0(ξ − ξ ')dξdξ '

C 0(ξ − ξ ')dξdξ '

(3.36)

(3.37)

Nella 3.36 il primo termine di sinistra è una costante ( C0(0) ), mentre il secondo termine è nullo poichè, essendo i due volumi v e vh disgiunti ed essendo a0 molto piccolo rispetto alla scala del lavoro, si ha che sempre ξ − ξ ' >> a0 e quindi C( ξ − ξ ' )= 0. Nella 3.37 il primo termine di sinistra è ancora la costante precedente, mentre il secondo termine, essendo i due volumi coincidenti, è un integrale non più nullo e può essere calcolato sulla base del ragionamento seguente. Poichè le dimensioni di v sono molto grandi rispetto ad a0 ( v>> a0 ) nell’integrale in questione quando ξ ‘ varia in v per una fissata posizione di ξ (fig. 3.51) si hanno contributi positivi soltanto in un intorno D( ξ ) molto piccolo di centro ξ e raggio a0.

per ξ '∈ D( ξ ) C 0 ( ξ − ξ ') ≠ 0

per ξ '∉ D( ξ ) C 0 ( ξ − ξ ') = 0

Fig. 3.51 - Integrazione di un effetto pepita.

68

La somma di questi contributi, trascurando gli effetti di bordo, è data da:

A=

v v

C 0(ξ − ξ ')dξdξ '= C 0(h )dh

L’ultimo integrale è esteso a tutto lo spazio e il valore di A è invariante rispetto a ξ . Quindi:

γ 0 v ( h) =

1 v2

v

Adξ =

A . v

A/v è dunque il contributo della struttura pepitica nella regolarizzazione. Sulla base di quanto è stato detto possiamo dare il contributo dell’effetto pepita nel calcolo della seguente espressione generale :

γ (v1, v 2) =

1 v1v 2

v1 v 2

γ (ξ − ξ ')dξdξ

che comprende, come caso parrticolare, i due termini di sinistra della 3.32. L’espressione precedente ricorre in molte operazioni geostatistiche. In essa i volumi v1 e v2 sono di dimensioni qualsiasi e localizzati in posizione qualsiasi. Il contributo della struttura pepitica è:

γ 0( v 1, v 2 ) = C 0( 0) − A

Mis ( v 1 ∩ v 2 ) . v 1v 2

(3.38)

SE v1 e v2 sono disgiunti la precedente vale: γ 0 ( v 1, v 2 ) = C 0 ( 0) . 1 Se v1 = v2 = v si ha: γ 0( v , v ) = C 0( 0) − A . v Se v1 ⊂ v2 si ha: γ 0 ( v 1, v 2 ) = C 0 ( 0) − A

1 . v2

3.7 Varianza di dispersione

In molte applicazioni, relative soprattutto a problemi di polluzione e di variabilità dei suoli, si è interessati al calcolo della dispersione spaziale di una variabile. Il suo valore generalmente dipende sia dalle dimensioni del supporto su cui essa è definita, sia dalle dimensioni del dominio entro cui si vuole studiare tale variabilità. In questo paragrafo vengono dati gli strumenti analitici per il calcolo della dispersione quando si conosce la funzione variogramma. Un classico problema è quando si conosce la dispersione di una variabile misurata su un supporto puntuale e si vuole calcolare il corrispondente valore per un supporto più grande.

69

Sia Z(x) una FA stazionaria di covarianza C(h) e variogramma γ ( h ) . Consideriamo un dominio di volume V all’interno del campo in cui la FA è definita (fig. 3.52). Per ogni punto x di V indichiamo con z(x) una realizzazione 1di Z(x) e con zv il valor medio di z(x) all’interno di V: zV =

1 V

V

z( x )dx.

1 ZV = V

Z (x)

V

Z ( x)dx

Fig. 3.52 - Definizione di un dominio V nel campo di definizione di Z(x).

La dispersione dei valori z(x) all’interno di V è data dalla varianza statistica: s2 =

1 V

V

( z( x) − z )2 dx. V

Il valore di s2 dipende dalla posizione di V nel campo di definizione di Z(x) e dai particolari valori di z(x). Se ora consideriamo in ogni punto x di V la variabile aleatoria Z(x) possiamo definire la variabile aleatoria: ZV =

1 V

V

Z ( x )dx.

(3.39)

Anche la dispersione S2 di Z(x) rispetto a ZV è una variabile aleatoria. Il suo valore atteso è chiamato varianza di dispersione di Z(x) rispetto in ZV ed è indicata con la notazione σ2 ( O / V ) , dove il simbolo O indiva che la dispersione si riferisce ad un elemento puntuale: 1 2 ( Z ( x ) − ZV ) dx . σ2 ( O / V ) = E S 2 = E (3.40) V V

1

Si ricorda che per realizzazione nel punto x della VA Z(x) si intende il particolare valore misurato o misurabile della variabile nel punto x.

70

In un modello stazionario σ2 ( O / V ) non dipende dalla posizione di V ma soltanto dalla forma e dimensioni di V ed è una caratteristica della FA Z(x). Sostituendo nella 3.40 l’espressione di Zv e sviluppando si ottiene: 1 V

σ 2 (O / V ) =

=

1 V

E [ Z 2 (ξ )] +

V

C (0) +

V

1 V2

1 V2

V V

V V

E [ Z (ξ ') ⋅ Z (ξ '')]dξ 'dξ ''−

C (ξ ' −ξ '')dξ 'dξ ''−

2 V

V

2 V

V

E [ Z (ξ ) ⋅ Z (ξ ')]dξ 'dξ

C (ξ − ξ ')dξ 'dξ

= C ( 0) + C (V ,V ) − 2C (V ,V ) = C ( 0) − C (V ,V ) = γ (V ,V )

(3.41)

L’espressione precedente, molto semplice da calcolare, da la varianza di dispersione di una variabile puntuale all’interno di un volume V. Essa dipende dalla funzione variogramma e dalla forma e dimensioni di V. Si noti (confronta l’espressione 3.24) che essa corrisponde alla quantità di cui si abbassa il sill passando da un variogramma puntuale ad uno regolarizzato di supporto V. Per V → ∞ si ha che: D2 ( O / ∞ ) = lim D2 ( O / V ) = lim γ (V ,V ) = C ( 0) . V →∞

V →∞

Quando C(0) esiste, cioè quando le ipotesi di stazionarietà sono verificate, il sill del variogramma, chiamato anche varianza a priori, coincide con la varianza di dispersione di un punto in un campo illimitato. Passiamo ora ad esaminare la dispersione all’interno di un volume V di una variabile non più puntuale ma definita su certo supporto. Suddividiamo il volume V in elementi di volume v. Se v<
1 n

n

Zvi =

i =1

1 n

n i =1

1 v

vi

Z ( x)dx =

1 V

V

Z ( x)dx.

Con procedura analoga a quella fatta in precedenza si può calcolare la varianza di dispersione di Zv all’interno di V:

σ 2 (v / V ) = E

1 n ( Zvi − ZV ) 2 . n i =1

71

Fig. 3.53 - Discretiz-zazione di V in una unione di elementi vi

Nella precedente sostituendo Zvi e Zv rispettivamente con: Zvi =

1 v

ZV =

1 V

vi

V

Z ( x )dx.

Z ( x )dx.

e sviluppando si ottiene: Nella precedente sostituendo Zvi e Zv rispettivamente con:

σ2 (v / V) =

−

1 n

n i =1

1 n

2 vV

n i =1

V V

1 v2

v v

E [ Z(ξ ' ) ⋅ Z(ξ ' ' )] dξ 'dξ ' '+

1 V2

V V

E [ Z(ξ ' ) ⋅ Z(ξ ' ' )] dξ 'dξ ' '+

E [ Z(ξ ' ) ⋅ Z(ξ ' ' ) ] dξ 'dξ ' '=

=

1 n {C (vi , vi ) + C (V ,V ) − 2C (vi,V )} n i =1

Essendo nella ultima espressione i primi due termini tra parentesi indipendenti dalla particolare posizione di vi, ed essendo: 1 n

n

C ( vi ,V ) = C (V ,V ) i =1

si ottiene:

σ2 ( v / V ) = C ( v , v ) − C (V ,V ) ed in funzione del variogramma: σ2 ( v / V ) = γ (V ,V ) − γ ( v , v )

(3.42) 72

Se v si riduce ad un punto, il termine γ ( v , v ) della espressione precedente si riduce a zero e quindi ritroviamo, come caso particolare, la 3.41. Se la funzione variogramma contiene una componente pepitica caratterizzata da : A = C (h)dh il suo contributo alla varianza di dispersione è dato, ricordando la 3.39 da:

C 0(0) −

A A 1 1 − C 0(0) − = A − V v v V

La linearità dell’espressione 3.42 conduce alla proprietà dell’additività della varianza di dispersione. Infatti siano v, V, G tre volumi tali che v ⊂ V ⊂ G , si ha:

σ2 ( v / V ) = γ (V ,V ) − γ ( v , v )

(3.43)

σ2 ( v / G ) = γ ( G , G ) − γ ( v , v )

(3.44)

σ2 (V / G ) = γ ( G , G ) − γ (V ,V )

(3.45)

Sostituendo nella 3.44 γ ( G , G ) con la sua espressione ricavata dalla 3.45, si ottiene:

σ2 ( v / G ) = σ 2 (V / G ) + γ (V ,V ) − γ ( v , v ) , cioè:

σ2 ( v / G ) = σ2 (V / G ) + σ2 ( v / V ) . Quindi la varianza di dispersione di una unità v nel dominio G è uguale alla somma della varianza di dispersione dell’unità v nell’unità inermedia V e della varianza di dispersione dell’unità V nel dominio G. Illustrazione di un esempio L’esempio che sarà illustrato è una elaborazione condotta su dati contenuti nella memoria: “Spatial Variability of soil Properties” di Goro Uehara et al. in Soil Spatial Variability, Proceedings of a Workshop of the ISSS and the SSSA, Las Vegas 1994. La mappa di fig. 3.54 rappresenta lo schema di campionatura di un sito dell’area di Kenana (Sudan). In essa sono stati prelevati 254 campioni si suolo su cui sono stati misurate diverse variabili tra cui l’ESP (Exchangeable Sodium Percentage), il cui 73

variogramma medio è riportato nella fig. 3.55a. L’aggiustamento più immediato è quello fatto con uno schema sferico di range di 3 Km ed un sill di 22. Si disponeva inoltre anche di 18 campioni di verifica ubicati in prossimità di punti campionati.

Fig. 3.54- Localiz-zazione della campionatura. a Kenana (da “Spatial Variability of soil Properties” di Goro Uehara et al. in Soil Spatial Variability, Proceedings of a Workshop of the ISSS and the SSSA, Las

Vegas 1994).

Ciò ha consentito una stima del variogramma anche a piccola scala. Il risultato è stato un variogramma con un sensibile effetto pepita (v. fig. 3.55b), per nulla ipotizzabile sul precedente grafico. Ciò dimostra quanto sia importante, quando si studia un fenomeno spaziale, campionare anche la piccola scala. Come si evince dalla figura, il nuovo variogramma è stato riaggiustato con uno schema sferico di range 4 e sill 16 ed un effetto pepita di 6.

Fig. 3.55 - Variogrammi sperimentali del ESP (da “Spatial Variability of soil Properties” di Goro Uehara et al. in Soil Spatial Variability, Proceedings of a Workshop of the ISSS and the SSSA, Las Vegas 1994).

74

Fig.3.56 - Andamento delle varianze di dispersione in funzione del supporto.

Utilizzando ognuno dei due variogrammi è stata calcolata la varianza di dispersione dell’ESP regolarizzato in funzione delle dimensioni del supporto, considerato come un quadrato di lato L. Nella fig. 3.56 sono riportati gli andamenti ottenuti con i due variogrammi. Come si nota la varianza diminuisce con l’aumentare di L. In particolare si può notare come il variogramma con effetto pepita (b) produce un forte calo della varianza di disper sione anche per piccoli supporti. Ciò è intuitivo: più è forte la variabilità a piccola scala, maggiore è l’effetto regolarizzante del supporto.

3.8 Il kriging ordinario

Una delle operazioni più comuni che vengono effettuate nell’ambito del trattamento dei dati spaziali, è la costruzione di carte tematiche, cioè carte georeferenziate relative a porzioni di aree geografiche, in cui è riportato, mediante un adeguato metodo di rappresentazione, l’andamento di una variabile di studio. La carta è normalmente costruita a partire dai valori della variabile misurati all’interno dell’area. Per esempio da una situazione di partenza espressa dalla fig. 3.57a, che rappresenta la posizione dei campioni ed i valori misurati di una variabile, si vuole costruire una carta ad isovalori come quella riportata nella fig. 3.57b. Si noti nell’esempio una distribuzione non uniforme della campionatura.

75

4955000

4950000

4945000

4940000

4935000

4930000

4925000 1640000

1650000

1660000

1670000

a)

b)

Fig. 3.57. Costruzione di una carta a partire dai punti di misura.

La carta ad isovalori, che nel gergo cartografico è anche chiamata carta vettoriale, è una delle tante forme di rappresentazione di una variabile geografica. Essa si ottiene non direttamente, ma successivamente ad una ricostruzione a maglia regolare della variabile tramite una operazione di stima (fig. 3.58a). Le linee di contorno si ottengono interpolando sui lati della maglia i valori da rappresentare (v. fig. 3.58b). Se ne deduce che la qualità di una carta è riconducibile alla qualità della stima che ha prodotto il grigliato di valori.

3.8.1 Gli stimatori lineari

Poichè i campioni sono raramente disposti a maglia regolare e, quand’anche lo fossero la densità della campionatura non sarebbe comunque sufficiente per operare un dettagliato tracciamento delle linee isovalori, la costruzione di un fitto grigliato di punti passa attraverso una operazione di stima, cioè una operazione mediante la quale viene attribuito in qualche maniera un valore ad una variabile in un punto in cui essa non è nota. La maniera con cui il valore viene calcolato definisce il tipo di stimatore.

76

Fig. 3.58 - Costruzione di una carta vettoriale: a) ricostruzione della variabile a maglia regolare; b) tracciamento per interpolazione delle linee di isovalori.

Questo tipo di operazione ha carattere locale in quanto la stima che si vuole effettuare non riguarda le caratteristiche generali (o caratteristiche globali) della variabile nel campo, ma riguarda appunto delle caratteristiche locali. Per tale ragione questo tipo di stima è chiamato nel gergo geostatistico stima locale. Gli stimatori più adatti e più usati per questo tipo di operazione sono quelli lineari. In essi il valore da attribuire ad un punto x0 del campo è calcolato mediante una combinazione lineare dei valori noti situati nelle vicinanze del punto da stimare, per es. entro un dominio circolare (vedi fig. 3.59). Tale dominio, che in realtà può avere qualsiasi forma e contenere alcune decine di punti (da 20 a 60) sarà chiamato d’ora in poi vicinaggio di stima.

77

Sia n il numero di punti noti che concorrono alla stima del valore in x0 e siano xα e z(xα) (α =1, n)1 rispettivamente le loro posizioni ed i corrispondenti valori della Fig. 3.59 Carattere dell’operazione di stima.

4955000

locale

4950000

variabile. Siano ancora z(x0) il valore sconosciuto in x0 e z* ( x 0 ) lo stimatore lineare considerato. Esso ha la seguente forma:

4945000

4940000

4935000

z* ( x 0 ) =

4930000

n

λαz ( x a )

(3.46)

α =1 4925000 1640000

Nella precedente λα sono i coefficienti della combinazione lineare, chiamati anche ponderatori dato che la combinazione lineare non è nient’altro che una media ponderata. 1650000

1660000

1670000

Proviamo ora ad usare gli strumenti geostatistici che abbiamo in possesso per caratterizzare questo stimatore e per definirne le condizioni ottimali in relazione anche agli stimatori lineari tradizionali. Poniamoci per semplicità nel quadro del modello stazionario, rimanendo inteso che le conclusioni alle quali arriveremo sono valide anche per i modelli quasi-stazionario e intrinseco. Sia pertanto Z(x) la FA stazionaria assunta per interpretare in senso probabilistico il fenomeno di studio e siano C(h) e γ(h) le funzioni covarianza e variogramma. Indicando con Z(x0) e Z(xα) le VA nei punti x0 e xα si ha per lo stimatore lineare: n

Z * ( x 0) =

λαZ ( x a )

(3.47)

α =1

A questa stima, come in genere ad ogni stima, è associato un errore, detto appunto errore di stima. Esso è dato dalla differenza tra il valore vero ed il valore stimato:

Z(x0) -

n

λαZ ( x a ) .

(3.48)

α =1

Esso continua ovviamente ad essere una VA. 1

La notazione xα ora introdotta sarà mantenuta anche nel seguito di queste dispense e starà sempre ad indicare la pozione dei punti di misura.

78

3.8.2 Correttezza della stima

Lo stimatore introdotto deve avere la fondamentale proprietà di essere corretto, cioè di essere di media nulla: E Z ( x 0) −

n

α =1

λαZ ( xa ) = 0

(3.49)

La precedente equazione diventa: E Z ( x 0) −

n

λαE Z ( x a ) = 0

(3.50)

α =1

In condizioni di stazionarietà E Z ( x ) = m e la (3.50) si riduce a:

m 1−

n

α =1

λα = 0 ,

che equivale alla seguente condizione sui ponderatori: n

λα = 1

α =1

3.8.3 Accuratezza della stima

La qualità della stima dipende dalla ampiezza degli errori (di stima). Essi sono caratterizzati da una legge di densità di probabilità che ha media nulla, se lo stimatore è corretto, e una dispersione che è responsabile dell’accuratezza della stima.(fig. 3.60). Più la funzione di densità è dispersa e più frequenti sono gli errori elevati. Come si può intuire, dal momento che si opera in un quadro probabilistico descritto da momenti del secondo ordine, non è possibile accedere a tale funzione, Fig. 3.60 - Distribuzione dell’errore di stima ma si è in grado, come si vedrà, di calcolarne la sua varianza, che viene chiamata varianza di stima ed è assunta quale grandezza per quantificare, in termini inversi, l’accuratezza della stima.

79

Useremo per tale varianza la notazione σ2s , dove il pedice s sta per stima. 2

σ = D Z ( x 0) − 2 s

2

= E Z 2 ( x 0) −

α

λαZ ( xα ) = E Z ( x 0) −

α

β

λαλβZ ( xα ) Z ( xβ ) − 2

λαλβC ( x a − x β ) − 2

= C ( 0) + α

β

α

α

λαZ ( xα )

λα Z ( xα ) Z ( x 0)

λαC ( x α − x 0 ) .

(3.51)

α

Sostituendo nella precedente C(h) con C ( 0) − γ ( h ) , poiché la somma dei ponderatori è uguale a uno, si ottiene la varianza di stima in funzione del variogramma:

σ2s = 2

λαγ ( x α − x 0) − α

λαλβγ ( x a − x β ) α

β

(3.52)

Si noti dall’espressione precedente che la varianza di stima non necessità dell’esistenza di C(0) e quindi della stazionarietà d’ordine 2. 3.8.4 Stimatori tradizionali

Esaminiamo dal punto di vista delle caratteristiche sopra specificate due degli stimatori tradizionali più usati: l’inverso delle distanze e l’inverso dei quadrati delle distanze. Inverso delle distanze: I ponderatori da attribuire ai campioni compresi entro il vicinaggio di stima hanno un peso proporzionale all’inverso delle distanze:

λα =

k , dα

con

1

k= i

1 di

Nella precedente di è la distanza di xi da x0

80

Inverso dei quadrati delle distanze I ponderatori da attribuire ai campioni hanno in questo caso un peso proporzionale all’inverso del quadrato delle distanze:

λα =

k , di2

con

1

k= i

1 di2

Come si può notare entrambi gli stimatori sono corretti, ed i valori dei ponderatori da attribuire ai campioni xα dipendono solo dalla loro distanza dal punto da stimare e non dalla loro reciproca posizione. Essi inoltre prescindono dalla variabilità spaziale del fenomeno di studio espressa dai variogrammi. Comunque, avendo in qualsiasi maniera calcolato dei ponderatori, se si dispone della funzione variogramma, è possibile calcolare mediante l’espressione (3.52) la precisione dello stimatore corrispondente. Diamo un esempio: si consideri lo schema di figura 3.61. Esso rappresenta la localizzazione di cinque punti di misura (1,2,3,4,5) e di un sesto, identificato da una crocetta, che è il punto in cui bisogna effettuare la stima. Considerate le distanze dei 5 punti da quello da stimare, i ponderatori derivanti dai due stimatori tradizionali sono riportati nella tabella 3.1. Nella ultima colonna della tabella è riportata la precisione della stima calcolata con la seguente funzione variogramma: γ (h) = 4.2 γ exp(h, 5) + 13.5γsph (h,18)

Si tratta del variogramma aggiustato per la variabile Cs137 campionata nell’area di Briansk, cittadina russa. Con questo variogramma la stima con l’inverso del quadrato delle distanze risulta essere leggermente più precisa. Fig. 3.61 - Configurazione di stima.

Tab. 3.1 Stimatore inverso distanze inverso quadrati distanze

λ1

0.1903 0.1747

λ2

0.1715 0.1418

81

λ3

0.1887 0.1718

λ4

0.2755 0.366

λ5

0.1739 0.1457

σs 2.813 2.829

3.8.5 Il Kriging Ordinario

Come si è visto, dato un set di ponderatori, si è in grado di calcolare, con l’ausilio del variogramma, la precisione della stima corrispondente. E’ immediato quindi porsi come problema quello di determinare quei ponderatori che danno luogo alla stima migliore, nel senso di una maggiore precisione. E’ un problema che si risolve minimizzando la varianza di stima, data dalla espressione (3.52), con il rispetto della condizione di correttezza. Il metodo di ottimizzazione è quello di Lagrange, che consiste nell’eguagliare a zero le n derivate parziali della (3.52) rispetto ai ponderatoti λα.:

∂σ 2s =0 ∂λ α

∀α = 1, n

con la condizione:

n

λα = 1.

α =1

Si ottiene:

∂σ s2 ∂ = 2 ∂λα ∂λα

λαγ ( xα − x 0) − a

α

β

λαλβγ ( xα − xβ ) + 2 µ (1 −

α

λα ) = 0 ,

∀α = 1,n

e svolgendo le derivate:

∂σ 2s = 2 γ ( x α − x 0) − 2 ∂λα

λβ γ ( x α − x β ) − 2 µ = 0 ,

∀α = 1,n

(3.53)

β

Nel seguito, quando ragioni di praticità lo richiederanno, adotteremo la notazione più compatta γαβ in luogo di γ ( x α − x β ) . L’espressione, ricordiamo, indica il valore della funzione variogramma per una distanza orientata che va, in questo caso, dalla posizione xα alla posizione xβ . Integrando la (3.53) con la condizione di correttezza si ottiene il seguente sistema:

λβγαβ + µ = γαo ,

∀α = 1,n

(3.54)

β

λβ = 1 β

Il sistema precedente è un sistema lineare di n+1 equazioni in n+1 incognite. Le incognite sono costituite dagli n ponderatori λα ed dal parametro µ, detto parametro di Lagrange. Lo stimatore associato è detto kriging ed il sistema (3.54) è chiamato sistema di kriging. Si dimostra che se le posizioni xα sono distinte, il sistema è sempre regolare, e quindi ammette sempre soluzione che è unica.

82

L’espressione (3.52) delle varianza di stima, che ormai possiamo chiamare varianza di kriging ed indicare con la notazione σk2, tenendo conto delle (3.54), può essere scritta nella forma più semplificata:

σ2k =

λβγβo + µ .

(3.55)

β

Scriviamo il sistema (3.54) in forma più esplicita:

λ1γ11 + λ2γ12 + ⋅ ⋅ ⋅ + λαγ1α + ⋅ ⋅ ⋅ ⋅ + λnγ1n + µ = γ10 λ1γ21 + λ2γ22 + ⋅ ⋅ ⋅ + λαγ2α + ⋅ ⋅ ⋅ ⋅ + λnγ2n + µ = γ20 λ1γ31 + λ2γ32 + ⋅ ⋅ ⋅ + λαγ3α + ⋅ ⋅ ⋅ ⋅ + λnγ3n + µ = γ30 ⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅ λ1γα1 + λ2γα + ⋅ ⋅ ⋅ + λαγαα + ⋅ ⋅ ⋅ ⋅ + λnγαn + µ = γα0 ⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅ ⋅⋅ λ1γn1 + λ2γn2 + ⋅ ⋅ ⋅ + λαγnα + ⋅ ⋅ ⋅ ⋅ + λnγnn + µ = γn0 λ1 + λ2 + ⋅ ⋅ ⋅ + λα + ⋅ ⋅ ⋅ ⋅ + λn = 1

(3.56)

Il sistema precedente, scritto in termini matriciali, diventa:

γ 11 γ 12 γ 21 γ 22

⋅⋅⋅

γ 1α γ 2α

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

γ αα

⋅⋅⋅

γαn

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

1 x λα = γα 0 ⋅⋅⋅ ⋅ ⋅⋅ ⋅ ⋅⋅

γ n1 γ n 2

⋅⋅⋅

γ nα

⋅⋅⋅

γnn

1

⋅⋅⋅

1

⋅⋅⋅

1

0

⋅⋅⋅

⋅⋅⋅

γα 1 γα 2 ⋅⋅⋅

1

1

⋅⋅⋅

1

⋅⋅⋅

γ 1n γ 2n

⋅⋅⋅

1

λ1 λ2

γ 10 γ 20

⋅⋅⋅

⋅ ⋅⋅

⋅ ⋅⋅

λn µ

(3.57)

γn 0 1

La matrice dei coefficienti non dipende dall’entità da stimare ma dipende esclusivamente dalla posizione reciproca dei punti di misura, cioè delle informazioni che si utilizzano per effettuare la stima, e dalla funzione variogramma. In particolare essa rappresenta la struttura spaziale dell’informazione. Il vettore dei termini noti è legato all’entità da stimare e descrive i rapporti spaziali tra questa ed i punti di misura. Il Kriging Raster Il sistema di kriging descritto sopra si riferisce alla stima di una variabile puntuale e perciò esso viene normalmente chiamato kriging puntuale. In molte applicazioni, soprattutto di carattere territoriale, si ha necessità di stimare le variabili su supporti areali, normalmente celle di forma quadrata o rettangolare. Tale necessità può sorgere per esempio nella integrazione delle stime in ambienti GIS, dove le elaborazioni delle carte avvengono su formati raster. Procedendo come per il caso puntuale possiamo dedurre le equazioni di kriging per la stima ottimale di una variabile su supporto raster.

83

Supponiamo che si debba stimare il valore della variabile Z estesa ad un volume v centrato in x0 a partire dal valore misurato della stessa su n punti di posizione xα. (fig. 3.62). L’entità da stimare è Zv(x0) = Zv0 =

1 v

v0

Z ( x )dx .

e lo stimatore lineare adottato Zv* ( x 0) è una combinazione lineare dei valori noti:

Zv* ( x 0) =

n

λαZ ( x a ) .

α =1

Fig. 3.62 Configurazione di una stima raster.

Si vuole una stima corretta e di varianza minima. Essendo E Z v ( x 0) = m , la condizione di correttezza si traduce, come per il caso della stima puntuale, alla condizione sui λα:

λα = 1. α

Per quando concerne l’altra condizione, si tratterà di minimizzare l’espressione della varianza di stima rispetto ai ponderatori tenendo conto della condizione di correttezza. La varianza di stima in questo caso è data da:

σ = D 2 s

2

1 v

v0

Z ( x )dx − α

λα Z ( xα ) = E

1 v

2 v0

Z ( x )dx − α

λα Z ( xα )

= 1 v2

v0 v0

=

1 v2

E [ Z ( x ) ⋅ Z ( x ')]dxdx '+

v0 v0

C ( x − x ')dxdx '+

α

α

β

β

λαλβ E [ Z ( xα ) ⋅ Z ( xβ )] − 2

λαλβC ( xα − xβ ) − 2

84

α

λα

α

1 v

λα

v0

1 v

v0

E [ Z ( x ) ⋅ Z ( xα )]dx

C ( x − xα )dx .

In funzione del variogramma essa si esprime:

σ2s = 2

λα

α

1 v

v0

γ ( x − xα )dx −

1 v2

v0 v0

γ ( x − x ')dxdx '−

α

β

λαλβγ ( xα − xβ ) .

Le equazioni di kriging si ottengono eguagliando a zero le derivate parziali della precedente tenendo conto della condizione di correttezza:

∂ σ s2 + 2 1 − ∂λα

α

λα

∀α = 1, n

=0

Si ottiene il seguente sistema:

β

λαγ ( xα − xβ ) + µ =

1 v

v0

∀α = 1,n

γ ( x − xa )dx

(3.58)

λβ = 1 β

con varianza di stima:

σ k2 =

β

λβ

1 v

v0

γ ( xβ − x )dx + µ −

1 v2

v0 v0

γ ( x − x ')dxdx '

Con scrittura più compatta si ha:

λαγαβ + µ = γα , vo

∀α = 1, n

β

(3.59)

λβ = 1 β

85

σ2k =

λβγβv + µ − γ ( v , v )

(3.60)

0

β

In termini matriciali il sistema si può scrivere:

γ 11 γ 12 γ 21 γ 22

⋅⋅⋅

γ 1α γ 2α

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

γαα

⋅⋅⋅

γαn

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

1 X λα = γ αvo ⋅⋅⋅ ⋅⋅⋅ ⋅⋅⋅

γ n1 γ n 2

⋅⋅⋅

γnα

⋅⋅⋅

γ nn

1

λn

γ nv

⋅⋅⋅

1

⋅⋅⋅

1

0

1

1

⋅⋅⋅

⋅⋅⋅

γα 1 γα 2 ⋅⋅⋅

1

1

⋅⋅⋅

1

⋅⋅⋅

γ 1n γ 2n

⋅⋅⋅

1

λ1 λ2

γ 1v γ 2v

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

o o

(3.61)

o

Come si può notare la matrice dei coefficienti è identica a quella della stima puntuale. Essa infatti dipende solo dalla configurazione dei punti di misura usati per la stima, che è identica nei due casi. E’ solo il termine noto che tiene conto della entità da stimare. Esso, per la generica linea α, esprime il valor medio che la funzione variogramma assume quando un estremo è fisso in nel punto di misura α e l’altro varia all’interno del volume da stimare v0 . Il suo calcolo normalmente viene effettuato numericamente discretizzando il volume in un insieme di punti e trasformando quindi l’integrale in sommatora: 1 v

v0

γ ( xα − x )dx =

1 N

N j =1

γ ( xα − xj )

Di solito una discretizzazione del volume in sedici punti a maglia regolare (v. fig. 3.63) è sufficiente a garantire una approssimazione accettabile. Fig. 3.63 Discretizzazione del volume da stimare

86

Prima di enunciare alcune importanti proprietà del sistema di kriging, esaminiamo come si comportano i ponderatori in un caso limite, quando il variogramma è un effetto di pepita puro, vale a dire un variogramma con range zero (o molto piccolo rispetto alla scala del lavoro):

γ ( h) =

0 se h = 0 C 0 se h ≠ 0

Consideriamo il caso di una stima puntuale e quindi facciamo riferimento al sistema 3.60. Consideriamo la generica riga α del sistema. Tutti i valori dei coefficienti valgono λβ = 1, la riga si riduce a: C0 eccetto γαα che vale 0. Pertanto, tenendo conto della β

C 0(1 − λα ) + µ = C 0 ,

ovvero:

λα =

µ C0

= C te

Quindi tutti i ponderati sono uguali e pari a 1/n. Il kriging si riduce ad essere una media aritmetica dei campioni selezionati per la stima. Si noti come l’assenza di correlazione spaziale rende i campioni tutti equivalenti ai fini della stima, qualunque sia la loro posizione nello spazio. Proprietà del Kriging. a) Il kriging è un interpolatore esatto Applichiamo le equazioni del kriging puntuale per la stima di un punto noto, supponendo che sia quello di posizione xβ.. Se si risolve il sistema di 3.57 con la regola di kramer si ha che: • quando l’incognita da determinare è λν ≠ λβ si ha che il determinante a numeratore dell’espressione risolutiva di Kramer vale zero, poiché la matrice corrispondente ha due colonne uguali (costituite dai termini: γαβ ( α = 1, n ) ): una già contenuta nella matrice dei coefficienti, l’altra costituita dalla colonna dei termini noti, che sostituisce la colonna corrispondente all’incognita. I ponderatori saranno sistematicamente nulli. • quando la incognita da determinare è λβ , i due termini che compaiono a numeratore e denominatore sono uguali e quindi il ponderatore vale 1. Il valore stimato perciò coinciderà con il valore vero: in questo senso il kriging è un interpolatore esatto.

87

b) Il kriging di una combinazione lineare coincide con la combinazione lineare del kriging dei suoi elementi. Analiticamente questa proprietà si esprime nella seguente forma:

[

]

*

p( x ) Z ( x )dx =

p( x ) Z * ( x )dx

(3.62)

Nella precedente si è utilizzata per maggiore generalità una misura di ponderazione p(x). L’asterisco indica la stima tramite kriging. Questa proprietà deriva dalla linearità delle equazioni di kriging. Diamo la dimostrazione facendo riferimento alla stima del valore medio di una variabile su una entità geometrica di volume v. Per la stima diretta del volume vale il sistema (3.61) dalla cui soluzione si ricavano i coefficienti, che per evitare ogni confusione, noteremo con λvα : Zv* =

1 Z ( x )dx v v

*

=

λα Z ( xα ) . (3.63) v

α

Se ora vogliamo stimare il valore di Zv come media delle stime effettuate sui punti interni di v dobbiamo considerare il sistema (3.57) applicato ad ognuno di essi. Sia x il generico punto all’interno di v (fig. 3.64). Notiamo con λα ( x ) i ponderatori corrispondenti. Fig. 3.64 - Stima dei punti interni al volume Integrando in v il sistema (3.57) si ottiene:

γ 11 γ 12 γ 21 γ 22

⋅⋅⋅

γ 1α γ 2α

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

γαα

⋅⋅⋅

γαn

1

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

γ n1 γ n 2

⋅⋅⋅

γnα

⋅⋅⋅

γ nn

1

⋅⋅⋅

1

λ n, v µv

γ n, v

⋅⋅⋅

1 0

⋅⋅⋅

⋅⋅⋅

γα 1 γα 2 ⋅⋅⋅

1

1

⋅⋅⋅

⋅⋅⋅

γ 1n γ 2n

1 1

λ 1, v λ 2, v

γ 1, v γ 2, v

⋅⋅⋅

⋅⋅⋅

⋅⋅⋅

X

λ α,v

=

dove si è posto:

λ α,v =

1 λα ( x )dx e v v

γ a, v =

1 γ ( xa − x )dx . v v

88

γα, v

1

Questo sistema coincide con il sistema (3.61): essi hanno la stessa matrice dei coefficienti e lo stesso termine noto. Avranno pertanto le stesse soluzioni:

λvα = λα , v . Questo implica l’uguaglianza tra le due entità:

Zv* =

α

λvα Z ( x α )

e 1 Z * ( x )dx = v v

α

1 λα ( x )dx ⋅ Z ( xα ) = v v

α

λ α , vZ ( xa ) .

La proprietà descritta in questo punto è una importante proprietà di congruenza del kriging, che gli stimatori tradizionali non posseggono.

c) Confronto con gli stimatori tradizionali Con riferimento alla configurazione (fig. 3.63) si possono calcolare e confrontare i ponderatori di kriging e la relativa varianza di stima con quelli degli stimatori tradizionali. La tabella che segue riprende la tab.1 mostrata in precedenza e la integra con i dati del kriging. Tab. 3.2 Stimatore inverso distanze inverso quadrati distanze Kriging

λ1

0.1903 0.1747 0.1981

λ2

0.1715 0.1418 0.2069

λ3

0.1887 0.1718 0.2212

λ4

λ5

0.2755 0.1739 0.366 0.1457 0.3297 0.441

σs 2.813 2.829 2.781

Si osservi l’avvenuta riduzione della varianza di stima. Un’altra forma di elaborazione è stata condotta per evidenziare la differenza tra gli stimatori tradizionali ed il kriging. Il caso è sempre quello del Cesio 137, da cui è stato tratto il variogramma utilizzato per il precedente confronto. Per ognuno dei tre stimatori sono stati stimati i valori del Cs137 nei 1978 punti dove il valore della variabile era noto, rimuovendolo di volta in volta. Un tale modo di procedere, che va sotto il nome di test kriging, è frequentemente usato in Geostatistica per il controllo del modello e dei parametri del modello. Relativamente ad ogni stima è stata elaborata la nuvola dei

89

correlazione valori veri/ valori stimati ed è stato costruito l’istogramma degli errori, calcolando la media e la varianza. Nella fig. 2.65 sono riportati questi risultati.

Z* 100

100

Z*

Z*

100

80

80

80

60

60

60

40

40

40

20

20

20

0

0 0

20

40

60

80

Z

f. rel.

0 0

100

20

40

60

80

Z

100

f. rel.

400

0

f. rel

400 300

300

200

200

200

100

100

100

0

0

-4

-3

-2

-1

0

1

2

3

4

5

a)

60

80

Z

100

0

-5

-4

-3

-2

-1

0

1

2

3

4

5

-5

-4

-3

-2

-1

0

1

2

Z-Z*

Z-Z*

m[Z-Z*]=0.0489 s2 [Z-Z*]=13.105 n=1978

40

400

300

-5

20

3

4

5

Z-Z*

m[Z-Z*]=0.060 s2[Z-Z*] =9.374 n=1978

m[Z-Z*]=0.017 s2[Z-Z*] =7.659 n=1978

b)

c)

Fig. 3.65 - Nuvola di correlazione veri/stimati ed istogramma degli errori: a) inverso distanze; b) inverso quadrati distanze; c) kriging.

I risultati ovviamente sono favorevoli al kriging: la varianza di stima è più piccola, l’istogramma degli errori è meno disperso e la retta di regressione, anche se non tracciata, si avvicina di più alla retta a 45°. Ancora un altro confronto con gli stimatori tradizionali: si osservi la fig. 3.66. In essa sono riportati cinque configurazioni di stima. La configurazione a) presenta quattro campioni localizzati sui vertici di un quadrato ed il punto da stimare occupa il centro del quadrato. La configurazione b) deriva dalla precedente rimuovendo il campione 4. Si passa poi dalla configurazione b) alla configurazione e) lasciando invariato il punto da stimare e avvicinado,via via fino a farli coincidere, il campione 1 al campione 2. Le stime sono state fatte impiegando un variogramma isotropo. Per la configurazione a) il kriging fornisce per i quattro pondaratori lo stesso valore, come è giusto che sia per simmetria. Anche i ponderatori tradizionali danno questo risultato. Per le restanti quattro configurazioni i metodi tradizionali attribuirebbero ai campioni uguali ponderatori. E ciò non sarebbe proprio. Nella configurazione b) il ponderatore 1 ed il ponderatore 3 devono essere uguali per simmetria, ma il ponderatore 2 è più piccolo degli altri due. Passando poi dalla configurazione a) alla configurazione d) il contributo dei ponderatori 1 e 2 , come si può intuire, deve diminuire, mentre quello del campione 3 deve aumentare. Come caso limite, quando il campione 1 ed il

90

campione 2 risultano coincidenti (configurazione d) il loro peso equivale a quello di un solo campione.

Fig. 3.66 - Configurazioni di stima con valore dei ponderatori di kriging.

Sensibilità ai parametri strutturali I risultati di un kriging consistono normalmente, per una data configurazione di stima, nel valore dei ponderatori e nella precisione della stima. Essi dipendono essenzialmente dal numero di punti di misura, dalla loro posizione geometrica, reciproca e rispetto alla entità da stimare, e dalla funzione variogramma. Per dare una idea della dipendenza dei risultati della stima dai parametri del variogramma abbiamo scelto una configurazione di stima a maglia regolare, che è riportata nella fig. 3.67. Essa si riferisce ad una stima puntuale in cui l’entità da stimare è posta al centro della configurazione e l’informazione è costituita da 12 punti di misura disposti e numerati come nella citata figura.

Fig. 3.67 - Configurazione scelta per la sensitività dei risultati del kriging ai parametri del variogramma.

91

Operando con un variogramma isotropo, le simmetrie sui campioni danno luogo a simmetrie sui ponderatori, per cui è comodo raccogliere i campioni in tre gruppi: (1, 2, 3, 4); (5, 6, 7, 8); (9, 10, 11, 12). Ai quattro campioni di ogni gruppo compete per simmetria un eguale valore del ponderatore, per cui sommando il loro contributo, e notandolo rispettivamente λ1, λ2, λ3, riduciamo a tre il numero dei ponderatori da considerare nell’analisi di sensitività.

Andamento dei ponderatori con il range L’analisi consiste nell’esprimere l’andamento dei pondaratori λ1, λ2, λ3 in funzione del range del variogramma. E stato considerato un variogramma di tipo sferico avente sill 1. e range a. Per svincolarsi dal particolare valore del range i ponderatori sono stati espressi in funzione del rapporto a/d, dove d è il lato della maglia elementare (v. la già citata fig. 3.67). Nella fig. 3.68 sono riportato i risultati dell’analisi. Si nota come per valori del range molto piccolo rispetto a d i tre ponderatori sono pressocchè uguali e valgono 1/3. Man mano che il valore di a aumenta (sempre rispetto a d) aumenta λ1

Fig. 3.68 - Andamento dei ponderatori con il range del variogramma

mentre λ2 e λ3 diminuiscono: aumenta in pratica il peso dei campioni più vicini al punto da stimare a scapito di quelli più lontani.

Andamento dei ponderatori con la frazione nugget del variogramma Con riferimento alla configurazione di fig. 3.68, si riporta nella fig. 3.69 l’andamento di l1, l2 e l3 in funzione del rapporto C0/(C0+C), dove C0 è la componente nugget del variogramma e C0+C il sill complessivo. Il calcolo è stato effettuato con un variogramma sferico di range pari a 5 volte il lato della maglia.

92

Fig. 3.69 - Andamento dei ponderatori con l' effetto pepita.

Relazione varianza di stima/densità di campionatura Sempre con riferimento allo schema di fig. 3.67 è stato calcolato il valore della varianza di stima in funzione del rapporto a/d. Il grafico corrispondente è riportato nella fig. 3.70. Si osservi come, partendo da bassi valori di a/d , all’aumentare di tale rapporto (cioè, per un fissato range, all’aumentare della densità di campionatura) si sconta subito un forte guadagno di precisione fino a valori di a/d di 3-4, divenendo poi il guadagno più modesto anche per forti aumenti della densità di campionatura. Anche in questo caso è stato considerato un sill unitario, per cui per ottenere il valore della varianza di stima quando il variogramma è uno schema sferico di sill c, basta moltiplicare per c il valore risultante dal grafico.

Fig. 3.70 - Andamento della varianza di stima in funzione della densità di campionatura.

93

3.9 Il kriging semplice Si consideri il caso una Funzione Aleatoria stazionaria Z(x) , in cui la media E[Z(x)] = m, costante in tutto il campo, è anche nota. Il kriging in queste condizioni è chiamato kriging semplice (KS) o kriging a media nota. Normalmente la media di una FA non si conosce, a meno di non disporre di innumerevole quantità di dati. Comunque il KS costituisce un importante riferimento teorico e metodologico e gode di importanti proprietà. Lavorare con una FA di media nota è lo stesso che lavorare con una FA di media nulla: infatti si può definire Y(x) = Z(x) – m che ha media nulla:

E[Y(x)] = E[Z(x)] - m = 0

(3.63)

Il kriging pertanto di Z(x0) a partire dalle informazioni Z(xα) (α =1,n) attraverso la stima di Y(x0):

si ottiene

Z * ( x0 ) = m + Y * ( x0 ) = m +

α

λα Y ( xα )

(3.64)

La condizione di correttezza per la stima di Y(x0) diventa:

E[Y ( x0 ) −

α

λα Y ( xα )] = 0

ovvero

E[Y ( x0 )] −

α

λα E[Y ( xα )] = 0

che per la (3.62), data la stazionarietà, è sempre verificata, qualsiasi sia l’insieme dei ponderatori {λ}. Per la varianza della stima vale la (3.51):

σ ks2 = C ( 0 ) +

λ α λ βC ( x a − x β ) − 2 α

β

λαC ( xα − x 0 ) α

che non è più equivalente alla (3.52) ottenuta dalla 3.51 sostituendo C(h) con C(0)-γ(h) e sfruttando la relazione

α

λα = 1

non più necessaria.

94

Derivando la (3.51) rispetto ai ponderatori ed eguagliando a zero si ottiene

β

λ β C ( x a − x β ) = C ( xα

−

(3.65)

∀ α = 1, n

x 0)

Esprimendo Y in funzione di Z la 3.64 diventa:

Z * ( x0 ) = m + Y * ( x0 ) = m +

α

λα [ Z ( xα ) − m] = m (1 −

Come si può osservare il complemento a 1 di α

α

λα ) +

α

λα Z ( xα )

λα costituisce il peso da attribuire alla

media. La varianza di stima espressa dalla 3.51 si può scrivere, tenendo conto delle equazioni 3.65, nella forma più esemplificata:

σ ks2 = C ( 0 ) −

λαC ( xα − x 0 ) α

Proprietà del kriging semplice: •

L’errore di stima Y(x0)-Y*(x0) è ortogonale sia con la generica variabile Y(xα) coinvolta nella stima che con la stima stessa Y*(x0). Infatti la covarianza:

E[(Y0 −

β

λ β Yβ )Yα ] = E (Y0Yα ) −

β

λ β E (Yβ Yα )

è nulla in virtù delle equazioni 3.65. Analogamente la covarianza:

E[(Y0 − = β

β

λβ Yβ )(

λβ C0 β −

α

α

β

λα Yα )] =

β

λβ E (Yβ Y0 ) −

α

β

λα λβ E (Yα Yβ ) =

λα λβ Cα β

è anch’essa nulla. Infatti, per ottenere l’eguaglianza dei due termini di destra, basta moltiplicare entrambi i termini della 3.65 per λα e sommare in α. Nelle precedenti per brevità sono state utilizzate delle notazioni più semplificate: Yα per Y(xα) e Cαβ per C(xα-xβ). •

Se la F.A ammette legge spaziale multigaussiana FY1 ,Y2 ,...,Yk ( y1 , y2 ,..., yk ) stazionaria (v.paragrafo 3.3) il kriging semplice della variabile Y(x0) calcolato a partire dalle variabili Y(xα) (α=1.n) e la corrispondente varianza di stima costituiscono la media e la varianza della legge gaussiana di Y(x0) noti Y(xα) (α=1.n) 95

•

−

E(Y0/y1,y2,…,yn) = Y*

−

2 Var (Y0/y1,y2,…,yn) = σ ks

La retta di regressione tra valori veri e valori stimati ha pendenza unitaria. Infatti, considerando la stima di Y nel punto x, la retta di regressione, che fornisce la migliore approssimazione lineare di E[Y(x)/Y*(x)], è data da:

Cov[Y * ( x), Y ( x)] * Yˆ ( x) = Y ( x) Var (Y * ( x))

Il coefficiente lineare è uguale a uno poiché Cov(Yx* , Yx ) =

λβ Cβx

β

e Var (Yx* ) =

α

β

λα λβ Cαβ =

β

λβ Cβx

in virtù della 3.65

Stima con una sola variabile Si supponga di voler effettuare la stima nel punto x0 a partire da una misura nel punto x0+h , il sistema 3.65, ridotto ad una sola equazione, diventa:

λh x C(0) = C(h) da cui λh = C(h)/C(0) = ρ(h) coefficiente di correlazione tra Z(x0) e Z(x0+h) . La varianza di stima σKS è data da C(0) – C(h)/C(0) x C(h) = C(0) [1-ρ2 (h)]. Quando h raggiunge il range, la varianza assume il valore massimo che vale C(0). Se nelle stesse condizioni la stima fosse stata effettuata con il kriging ordinario, il ponderatore sarebbe stato 1 e la varianza di stima σKO sarebbe stata 2γ(h) = 2[C(0)C(h)], sistematicamente, eccetto che per h=0, più alta di σKS. Il rapporto σKS/σKO risulta essere infatti:

σ KS 1 C ( h) = 1+ σ KO 2 C (0 ) Il rapporto assume valore uno per h=0 e diminuisce progressivamente fino a 0.5 quando h raggiunge il range. Il KO, che è operato a media non nota, implica la stima della media e pertanto la varianza di stima è più alta di quella del KS dovendo tener conto dell’incertezza sulla media.

96