Mappatura del Radon in Alto Adige:
un’analisi di tipo geostatistico
Luca Verdi
Agenzia Provinciale per l’Ambiente di Bolzano, via Amba Alagi 5, 39100 Bolzano,
[email protected]
Stefano Pegoretti
Dipartimento di Fisica, Università degli Studi di Trento, via Sommarive 14, 38050 Povo (TN),
[email protected]
Riassunto
Vengono presentati i risultati ottenuti da una prima applicazione dei metodi geostatistici all’analisi della
distribuzione spaziale di radon relativamente al territorio altoatesino. Il dataset di partenza è costituito da oltre
4000 valori di concentrazione; analizzando i fattori descrittivi che accompagnano ogni misura, si è reso tale
dataset più omogeneo, riducendone al contempo la numerosità. La struttura spaziale dei dati è stata investigata
mediante analisi variografica, rivelando la presenza di un elevato “effetto nugget”; la “cross-validation” condotta
sui risultati ottenuti con l’approccio del Kriging Ordinario ha messo in luce un forte effetto di smoothing legato alla
procedura di interpolazione stessa. Al fine di considerare il più possibile l’intera variabilità del fenomeno, sono
state condotte delle Simulazioni Gaussiane Sequenziali: questo tipo di simulazioni consente di ottenere delle
Mappe di Probabilità di eccedere un determinato valore di soglia, con la possibilità di definire regioni spaziali a
elevata probabilità di alte concentrazione di attività di radon, come previsto dalla legge italiana.
A) INTRODUZIONE
“La geostatistica offre un insieme di strumenti deterministici e statistici volti alla comprensione e
modellizzazione della variabilità spaziale” {1}; così, sarebbe opportuno ricorrere a tale disciplina in
tutte quelle situazioni per le quali la continuità spaziale del fenomeno in esame si ritiene sia l’aspetto
principale, o quantomeno uno dei principali.
La letteratura si riferisce spesso alla geostatistica come alla teoria delle “variabili regionalizzate”,
ovvero dati che risultano spazialmente distribuiti e spazialmente strutturati. Fu Matheron, nel 1963, a
introdurre per primo questo termine, per mettere in luce i due aspetti principali che caratterizzano
questo tipo di dati: da un lato, l’aspetto “casuale”, che rende conto della variabilità spaziale su piccole
scale e di eventuali irregolarità locali; dall’altro, un aspetto più “strutturato”, che rende conto della
variabilità su grande scala e può evidenziare la presenza di tendenze, come trend e/o particolari
direzionalità del fenomeno.
Rispetto all’approccio statistico tradizionale, la geostatistica offre la possibilità di trattare allo stesso
tempo entrambi questi aspetti contrastanti; inoltre fornisce gli strumenti adeguati per gestire le
caratteristiche intrinseche di questo tipo di dati, che ne possono complicare in maniera sensibile la
trattazione:
ƒ generalmente si ha a disposizione una realizzazione unica del fenomeno per ogni
localizzazione campionata;
ƒ i dati spaziali non possono essere considerati indipendenti tra loro ed estratti dalla stessa
distribuzione di probabilità.
La geostatistica tradizionale si fonda sulla correlazione spaziale basata su una statistica “a due
punti” {2,3}; tuttavia recentemente si stanno sviluppando approcci che considerano una statistica a più
di due punti e approcci legati alla teoria bayesiana della massima entropia {4}.
1) LA LOGICA DELLO STUDIO GEOSTATISTICO
Prima di tutto risulta necessario disporre di uno strumento adatto a riconoscere, descrivere
e―possibilmente―modellare la variabilità dei dati; tra tutte le misure analitiche della
correlazione/continuità spaziale disponibili al giorno d’oggi, quella tradizionale e storicamente più
utilizzata è il variogramma γ(h), che dal punto di vista fisico descrive, al variare della distanza h,
quanto due misure sono tra loro scorrelate.
Una volta in possesso di un modello analitico per la correlazione spaziale, questo verrà
implementato in un secondo strumento atto alla stima/predizione della variabile di interesse nelle
localizzazioni che non sono state campionate. Limitandosi comunque all’insieme dei modelli lineari1, la
geostatistica ha sviluppato un’intera famiglia di stimatori, chiamati in gergo “Kriging”, che
appartengono ai cosiddetti stimatori BLUE(P) ― Best Linear Unbiased Estimators/Predictors, dove
Best indica che lo stimatore è scelto in modo da minimizzare la varianza della stima, Linear che si
tratta di modelli lineari, Unbiased che, nel senso del “valore di aspettazione”, non è affetto da errori
sistematici.
I pesi che la procedura di kriging assegna alle misure impiegate nella stima sono determinati al fine
di considerare in modo esplicito la distribuzione spaziale locale delle misure stesse attorno al punto di
stima; sono inoltre influenzati dal modello analitico di variogramma scelto nella prima fase dello studio
sulla base dell’analisi dei campionamenti disponibili.
B) IL DATASET IMPIEGATO
1) DESCRIZIONE DEL DATASET COMPLETO
Il dataset completo messo a disposizione dall’APPA di Bolzano è costituito da 4064 misure di
concentrazione di attività di radon indoor georeferenziate, per ognuna delle quali sono state raccolte
numerose informazioni, sia di tipo nominale che numerico, al fine di caratterizzare i principali fattori
che potrebbero influenzare la misura stessa, come la geologia della zona, il tipo e la qualità della
costruzione sede della misura, il tipo di locale, la stagione di esposizione, ecc.
Con lo scopo di rendere più affidabile il valore di concentrazione ottenuto, i valori ricavati da
dosimetri del tipo LR-115 (93% del totale) sono stati opportunamente corretti per un fenomeno noto
come “saturazione” {7,8}: questo tipo di intervento risulta particolarmente significativo nei casi di
elevati valori di concentrazione, conducendo così a significativi cambiamenti nella coda superiore
della distribuzione di probabilità che caratterizza i dati e che, come riconosciuto anche in letteratura, è
stata assunta di tipo log-normale [cfr. (fig. 1)].
Figura 1
(b)
(a)
Istogramma (a) e probability-probability plot (b) per i valori di concentrazione di radon indoor corretti per il
fenomeno della saturazione; la funzione di probabilità cumulativa teorica di riferimento è di tipo log-normale.
2) COSTRUZIONE DEL DATASET PER LO STUDIO GEOSTATISTICO
In linea di principio l’approccio di tipo geostatistico richiede che il dataset sia il più omogeneo
possibile, ovvero che tutte le misure appartengano alla stessa popolazione di riferimento. Basandosi
su precedenti esperienze e su informazioni di tipo fisico, per questo studio preliminare si sono
selezionate solo le misure condotte nel semestre invernale al piano terra; questa restrizione sui dati
ne ha ridotto la numerosità a 2578 (63% del totale).
Inoltre si è scelto di applicare una procedura di jackknife su tali dati, ovvero di dividere l’intero
dominio di studio in celle di dimensione (12 x 7.7) km2, per ognuna delle quali è stato estratto a caso
un valore; in questo modo, si è ottenuto un secondo dataset (costituito da 85 campioni) con la stessa
copertura spaziale di quello principale. Il dataset così prodotto è stato impiegato successivamente per
la validazione del modello geostatistico proposto, basato e costruito riferendosi unicamente al primo
1
Modelli per i quali un valore non noto della variabile viene stimato mediante una combinazione
lineare (eventualmente caratterizzata da opportuni pesi) di un fissato numero di dati misurati.
Figura 2
-3
Estimated Rn Concentration (OK) [Bq m ]
4000
3500
3000
2500
2000
A
1500
1000
500
B
0
0
500
1000
1500
2000
2500
3000
3500
4000
-3
Observed Rn concentration [Bq m ]
(a)
166850
168750
3907 (4)
166800
74 (4)
2474 (3)
326 (6)
1232 (1)
point A [135]
166650
808 (8)
166600
409 (7)
168250
230 (3)
81 (9)
1530 (9)
166500
1212 (10)
1706 (1)
168000
point B [3794]
50 (8)
166550
194 (6)
473 (5)
y coord [m]
y coord [m]
168500
1536 (2)
166750
166700
59 (7)
3618 (5)
319 (2)
167750
69 (10)
166450
611700
611750
611800
611850
x coord [m]
611900
611950
612000
167500
695000
695500
696000
696500
697000
x coord [m]
(b)
(c)
(a) risultati di uno studio di cross-validation; (b) e (c) zoom locale dei valori misurati di concentrazione e
relativa distribuzione spaziale dei campionamenti attorno ai due punti “patologici” A e B; per ogni punto, si
-3
riportano la concentrazione espressa in Bq m e, tra parentesi quadre, il valore reale del punto “patologico”.
La concentrazione stimata per tale punto si può ricavare da (a).
dataste, al quale vengono sottratti gli 85 valori necessari per il jackknife (costituito quindi da 2493
misure).
C) MAPPATURA MEDIANTE KRIGING ORDINARIO (OK)
Inizialmente per indagare sulla struttura spaziale dei dati, sono stati calcolati alcuni variogrammi
omnidirezionali mediante le librerie GSLib {1}; successivamente, riconosciuta tale struttura, si è ricorsi
al software Geostat Office {4} al fine di mettere in luce eventuali anisotropie spaziali. In relazione alla
rumorosità dei dati, si è dovuto ricorrere a una misura della correlazione spaziale non convenzionale,
il cosiddetto Pairwise Relative Variogram {1,3}, impiegato per la costruzione della superficie
variografica2 sull’intera scala regionale (circa un centinaio di km). Questo tipo di analisi ha rivelato la
presenza di anisotropia sia geometrica che zonale, ed entrambe sono state inglobate nel modello di
variogramma costruito e scelto in base a studi di cross-validation e jackknife.
Al fine di prendere in considerazione la non-stazionarietà del fenomeno (presenza di
eteroschedasticità―ovvero la varianza locale dei dati non è costante sull’intero dominio e di effetto
proporzionale―ovvero media e varianza locali sono linearmente dipendenti) le mappe per le stime del
valor medio di concentrazione di radon sono state ottenute ricorrendo all’algoritmo del Kriging
Ordinario con un’area di ricerca dei punti di “vicinaggio” mobile: solo i campionamenti che ricadono in
un cerchio di raggio prefissato centrato sul punto di stima vengono impiegati nel calcolo dei pesi del
kriging e conseguentemente nella determinazione della stima di media e varianza. In questo modo, si
ha la possibilità di recuperare una sorta di stazionarietà locale.
2
Un grafico polare bidimensionale in cui vengono assemblati diversi variogrammi direzionali: in questo
modo, si ha un’idea “visiva” di come la struttura spaziale eventualmente cambi con la direzione.
Questa e altri tipi di analisi hanno messo in luce la presenza di un elevato “effetto nugget”, ovvero
elevati valori della varianza su piccole scale; questo fenomeno riduce inevitabilmente l’efficacia di
qualsiasi tipo di kriging cui si intenda ricorrere. Come esempio, riportiamo in (fig. 2b) e (fig. 2c) due
tipiche situazioni “patologiche” riscontrabili nel dataset utilizzato che mettono in luce nella pratica
questa “indesiderabile” variabilità su piccola scala.
In (fig. 2a) sono riportati i risultati di uno studio di cross-validation con lo scopo di mettere in luce
l’effetto di smoothing legato all’interpolazione con OK: i valori elevati di concentrazione vengono
sistematicamente sottostimati, mentre quelli più bassi sono sovrastimanti; la pendenza della retta di
interpolazione (ottenuta mediante minimi quadrati) è pari a 0.24. Questo effetto è tuttavia comune a
qualsiasi tipo di interpolatore e in linea di principio non può essere evitato.
Figura 3
-3
Mappa ottenuta mediante OK per il valore medio della concentrazione di radon indoor (valori espressi in Bq m );
la griglia di interpolazione ha una risoluzione pari a 1 km in entrambe le direzioni; il raggio che definisce l’area di
ricerca dei punti di vicinaggio è pari a 5 km e il numero massimo di punti coinvolti nella stima è pari a 5; i punti
sulla mappa rappresentano le localizzazioni campionate.
La fig. 3 presenta una mappa della distribuzione spaziale della stima del valor medio di
concentrazione di radon indoor: una mappa di questo tipo rappresenta uno strumento atto alla
“descrizione” delle proprietà spaziali medie del fenomeno piuttosto che alla “previsione” accurata del
fenomeno.
D) MAPPATURA MEDIANTE SIMULAZIONI STOCASTICHE (C-sGs)
Nella pratica spesso si ha l’esigenza di disporre di uno strumento mediante il quale poter dare delle
risposte a domande come: Qual è la probabilità di superare un certo valore di soglia in una
localizzazione non campionata? Tra le varie possibilità offerte dalla geostatistica, la scelta è ricaduta
sull’approccio delle Simulazioni Stocastiche Gaussiane Sequenziali (sGs) condizionate ai
campionamenti; ogni simulazione deve rispettare i valori misurati nelle localizzazioni campionate.
L’idea che sta alla base di questa procedura è quella di sviluppare una sorta di metodo Monte
Carlo spaziale in grado di generare realizzazioni multiple ed equiprobabili del fenomeno; ognuna di
queste realizzazione viene indicata come immagine stocastica. Inoltre le simulazioni condizionate
possono essere usate come un modello realistico della variabilità spaziale {4}.
Come primo tentativo si è ricorsi all’approccio delle simulazioni gaussiane in quanto nella pratica il
modello gaussiano multivariato dovrebbe essere la scelta primaria per la simulazione di variabili
continue {1}. Sottolineamo infine come le simulazioni stocastiche abbiamo, in linea di principio, il
vantaggio di poter riprodurre l’intera variabilità del fenomeno sotto esame.
Come richiesto dall’algoritmo delle sGs, l’analisi variografica e le successive simulazioni sono state
condotte su dati NScored: i dati “grezzi” sono stati trasformati, secondo una funzione densità di
probabilità cumulativa (cdf) empirica ricavata dai dati stessi, in modo da ottenere una distribuzione
normale standard (media nulla e varianza unitaria). L’assunzione di una distribuzione normale
bivariata è stata verificata e ai risultati ottenuti in fase di simulazione è stata successivamente
applicata una trasformazione inversa secondo la stessa cdf empirica usata per la trasformazione
diretta.
In fase di simulazione si è deciso di ricorrere a OK con area di ricerca dei punti di vicinaggio mobile
per ottenere la stima della media locale della funzione di distribuzione normale cumulativa, mentre per
la stima della varianza locale della stessa distribuzione si è applicato l’algoritmo del Simple Kriging
(come raccomandato dalla letteratura). Al fine di evitare possibili artefatti legati ai valori casuali estratti
dalla distribuzione normale, si è adottato un approccio di tipo “multigrid”. L’intera parte computazionale
è stata eseguita con Geostat Office.
Rispetto alle mappe prodotte con tecniche di kriging o altri interpolatori, il reale vantaggio delle sGs
risiede nella possibilità di produrre mappe di probabilità [cfr. (fig. 4)]: dopo aver elaborato un numero
sufficiente di immagini stocastiche, si ottiene uno strumento affidabile e adeguato per identificare le
aree che manifestano un’elevata probabilità di superare un prefissato valore di soglia―come richiesto,
ad esempio, dalla legge italiana {5} e da direttive a livello comunitario {6}.
Figura 4
-3
Probabilità di superare una concentrazione di radon indoor pari a 400 Bq m ; la risoluzione della griglia di
simulazione è pari a 1.5 km in entrambe le direzioni; il raggio che definisce l’area di ricerca dei punti di vicinaggio
è pari a 10 km e il numero massimo di punti coinvolti nella stima è pari a 5 per i valori reali e a 3 per quelli
simulati; la mappa riprodotta è il risultato dell’elaborazione di 250 immagini stocastiche.
E) CONCLUSIONI
Nel presente lavoro si sono applicati gli strumenti messi a disposizione dalla geostatistica al fenomeno
“radon indoor” che, almeno nel caso esaminato, si caratterizza per (i) un’elevata variabilità dei dati su
scale spaziali differenti e per (ii) una complessità orografica del territorio di studio; entrambi questi
fattori hanno delle notevoli e inevitabili ripercussioni sulla qualità delle stime che si possono ottenere.
Tra tutti gli strumenti disponibili al giorno d’oggi, si sono scelti inizialmente due differenti approcci di
tipo geostatistico, basando la scelta sulla relativa semplicità (anche teorica) degli strumenti, sulla loro
affidabilità e sull’obiettivo dell’analisi; in particolare:
ƒ il metodo del Kriging Ordinario con area di ricerca dei punti di vicinaggio mobile si è
rivelato come un buono strumento che consente di trattare fenomeni non-stazionari e in
grado di produrre delle mappe delle proprietà medie del fenomeno sotto esame sull’intera
scala di interesse (un centinaio di km);
ƒ le Simulazioni Stocastiche Gaussiane e Sequenziali consentono teoricamente di prendere
in considerazione l’intera variabilità manifestata dai dati, evitando così l’effetto di
smoothing comune a tutti gli algoritmi di interpolazione (kriging incluso); elaborando un
numero sufficiente di immagini stocastiche, è possibile ottenere mappe di probabilità di
superare un fissato valore di soglia per la concentrazione di radon indoor, e
conseguentemente di identificare sul territorio delle possibili “radon prone areas”.
Infine, anche se i dati relativi a misure di radon possono risultare―e spesso risultano―molto
complessi per la loro stessa natura, questo lavoro dimostra che l’approccio geostatistico può essere
applicato anche a questo tipo di dati, eventualmente ricorrendo a misure della correlazione spaziale e
a strumenti non-convenzionali, quali ad esempio variogrammi relativi.
Come passo successivo, nell’ottica di migliorare la qualità delle stime ottenute, si procederà a
modificare il dataset in questo modo: misure appartenenti a un medesimo edificio (ma, per esempio, a
differenti stanze) saranno aggregate in un singolo valore di concentrazione3. La procedura proposta
dovrebbe, in linea di principio, ridurre l’effetto nugget e probabilmente l’intera varianza, mettendo di
conseguenza a disposizione un dataset più “leggibile”. A questo nuovo dataset verranno applicate
anche tecniche statistiche e geostatistiche più avanzate.
Bibliografia
{1} Deutsch, C. V. and Journell, A. C. GSLIB. Geostatistical Software Library and User’s Guide.
Applied Geostatistics Series (Oxford University Press), Second Edition 1998.
{2} Isaaks, E. H. and Srivastava R. M. An Introduction to Applied Geostatistics. Applied Geostatistics
Series (Oxford University Press), 1989.
{3} Goovaerts, P. Geostatistics for Natural Resources Evaluation. Applied Geostatistics Series (Oxford
University Press), 1978.
{4} Kanevski, M. and Maignan M. Analysis and Modelling of Spatial Environmental Data. (EPLF
Press), 2004.
{5} Decreto Legislativo 26-5-2000, n. 241, Attuazione della direttiva 96/29/Euratom in material di
protezione sanitaria della popolazione e dei lavoratori contro i rischi derivanti dalle radiazioni
ionizzanti.
{6} European Commission, Radiation Protection 88, Recommandations for the implementation of Title
VII of the European Basic Safety Standards Directive (BSS) concernine significant increase in
exposure due to natural radiation sources, 1997.
{7} Campi, F., Caresana, M., Ferrarini, M., Garlati, L., Palermo, M. and Rusconi R. Uncertainty
evaluation of Radon Measurements with LR115 Detector and Spark Counter. Rad. Prot. Dosim.,
Vol III, No. 1, pp. 59-54, 2004.
{8} Pegoretti, S. Mappatura del Radon in Alto Adige: un Approccio Geostatistico. Tesi di Laurea
Specialistica, Università degli Studi di Trento, 2005.
3
Attualmente, ci sono infatti molte situazioni per le quali campionamenti molto vicini tra loro
manifestano valori di concentrazione molto differenti.