Konstantin Krivoruchko
Environmental Systems Research Institute, Inc.
380 New York St., Redlands, CA 923738100, USA
kkrivoruchko@esri.com
Abstract
A large amount of data concerning the ecological state of Belarus after the Chernobyl accident was analyzed using geostatistics. This data included soil contamination by longlived radioisotopes ^{137}Cs, ^{90}Sr, ^{238240}Pu, and ^{241}Am; dose load estimates on the population during the first days after the Chernobyl accident; ^{137}Cs food contamination and estimation of the internal dose in 1993; and thyroid cancer morbidity among adults and children based on data collected from 1986 to 1995.
Currently, in order to integrate powerful methods of spatial data interpolation into a geographic information system (GIS), users are required to use statistical packages to process and then save the results of interpolation into a GISsupported format. To analyze complicated environmental problems and to present the results of an analysis on a map requires a group of specialists who are familiar with both GIS techniques and complicated statistical software. In this article we discuss an approach of combining advanced spatial statistics and GIS. This combination was first presented in the GIS MapStudio (Krivoruchko et al., 1997; Maignan, Krivoruchko, 1997) software package, which integrated both modern visualization techniques and geostatistical methods of spatial data analysis. These methods are currently being considered for further development and implementation into ESRI’s GIS software line.
This paper highlights the advantages of using geostatistical methods for processing environmental data. The following geostatistical approaches are used for mapping: simple, ordinary, indicator, probability, and disjunctive krigings, as well as ordinary and indicator kriging modifications for binomial data. The applicability of the different data processing methods is discussed. The data and the visualization techniques presented in this paper can help reveal powerful patterns for decision making and policy planning.
Introduction
The Chernobyl accident affected Belarus more than any other country. Ironically, Belarus does not have a single nuclear power plant; however, a ring of nuclear power stations surrounds it including the Ignaline station in Lithuania, the Smolensk station in Russia, and the Chernobyl station in Ukraine. The windroses on April 2630, 1986, was such that about 70 percent of the radioactive dust from Chernobyl fell on Belarus. The analysis and visualization of environmental and epidemiological data, therefore, are critical for predicting public health issues. Many interpretations of the effects from the Chernobyl accident appeared in postChernobyl years. The effects of radiation 12 years following the accident, however, are yet to be fully explored and are therefore poorly understood. The situation has been exacerbated by a lack of experience using modern statistical interpolation methods together with a GIS. Such techniques can reveal interrelationship between spatial distributions of radioactive contamination and disease incidence rates.
The Chernobyl accident released about 1.85× 10^{18} Becquerel of radioactive material. Each Becquerel represents one radioactive decay event per second. In 1958, Andrey Sakharov described "nonthreshold effects," by which every radioactive particle released had a statistical probability of doing damage either to the DNA of a cell or to the immune system through lowlevel internal radiation caused by ingesting such particles. He also predicted that radiation would accelerate the mutation of microorganisms, leading to the inference that persons with damaged immune systems would in time succumb more easily to these new strains (Sakharov, 1958). The epidemiological studies carried out in Belarus have revealed an increase in the incidence of thyroid and lung cancers, cardiovascular pathologies, and pregnancy anemias and an increase in mortality due to birth defects occurring in the contaminated regions (Lutsko et al., 1996).
During the first days after the accident, residents of Belarus absorbed the majority of their radiation dose through their thyroid glands via inhalation of contaminated air and, more importantly, consumption of contaminated foodstuffs, mainly cow's milk. Levels of contamination by shortlived radionuclides, in particular ^{131}I, were so high that the corresponding exposure of millions of people was qualified as "iodine shock". This will result in longterm health problems for the population. The epidemic of childhood thyroid cancer is the first indisputable long term health aftereffect of the accident.
At present (and over the next decades), the main hazard comes from ^{137}Cs (^{90}Sr, ^{238240}Pu, and ^{241}Am also play a significant role at short distances from the reactor). Today, internal exposure from intake of ^{137}Cs contaminated food contributes to more than half of the whole radiation dose received by the population. Income of the inhabitants of villages in Belarus does not afford them access to "clean" (non contaminated) food. They consume vegetables, potatoes, and milk produced on their own contaminated personal properties as well as mushrooms and berries from nearby forests.
The main goal of investigating the Chernobyl accident is to estimate the risk of health hazard to the population, which was irradiated and continues to live in the contaminated territory, based upon information about the total exposure of the population. In the current article we attempt to:
Geostatistics are statistical methods used to describe spatial relationships among sample data and to apply this analysis to the prediction of spatial and temporal phenomena. They are used to explain spatial patterns and to interpolate values at unsampled locations. Geostatistics have traditionally been used in the sphere of geosciences: meteorology, mining, soil science, forestry, fisheries, remote sensing, and cartography. Geostatistical techniques were originally developed by Soviet scientists for meteorological data predictions. The first book with complete explanations about simple and ordinary kriging and cokriging techniques was published in Leningrad in 1963 (Gandin, 1963). According to this book, the original name of the technique is objective analysis. Geostatistical techniques were later successfully applied to mining and other disciplines.
In addition to the description of spatial patterns and interpolation of data, important components of geostatistical analysis include the integration of secondary data in prediction algorithms through the use of cokriging and error analysis. The cokriging technique allow one to compute an optimum estimation while considering the relationships between several associated variables. Using error analysis, validation, and crossvalidation, it is possible to assess the performance of the interpolation as well as to reveal errors in the source data.
Sampling and mapping in the earth sciences are complicated by complex spatial and temporal variations. The structure and intensity of patterns being sampled often cannot be determined or predicted reliably with deterministic models because of uncertainties in the data and the phenomena under investigation. The best we can do using interpolation and estimation methods is to be as objective as possible and to consider the interrelationship of the data under investigation.
Deterministic approaches to interpolation (e.g. trend surface, inverse distance weighting, triangulation, and splining) are based upon a priori mathematical models of spatial variation and can lead to smooth but inaccurate maps. This is because error is an inherent part of the sampling process. In practice, error can not be eliminated but only minimized. Therefore, in most cases one cannot produce a representative map of estimated values in unsampled locations with these techniques.
Geostatistical processes are comprised of the three major components of regionalized variable analysis: variographical analysis, prediction making, and error analysis.
During structural analysis, spatial (auto)correlation can be analyzed using covariance and semivariogram. After structural analysis, predictions at unsampled locations are made using kriging (i.e. transposition of multiple linear regression into a spatial context).
The value of the error variation shows the uncertainty of the data at the considered point. Kriging variance is higher in poorly informed zones. Simple and ordinary kriging variance has a limited usefulness in the case of nonGaussian distribution. In this case the map of error of estimation associates with the data density and kriging variance equals to average variance of the selected neighborhood. If normal score transform had used and the bivariate normality was adopted, we can correctly calculate the confidence interval of estimation based on the features of the normal distribution. Thus, map of error estimation will provide an accurate measure of the local estimation precision.
This fact makes it possible to obtain maps of the probability of exceeding some predicted level by the variable under investigation (i.e., kriging provides the user with a tool that can help in the decision making process) as well as to draw the isolines with fixed significance.
For producing riskqualified maps, the nonparametric geostatistical algorithms such as indicator kriging, indicator cokriging, and probability kriging are primarily recommended. They were developed to process data of highly variant phenomena without having to trim off important highvalued data.
Validation allows one to randomly select parts of the database, perform kriging, and then compare results of estimation with true values. One more procedure of error analysis, in addition to validation and mapping errors of estimation or errors of indicators, is crossvalidation. This method consists of analyzing and displaying estimated values calculated under the assumption that particular sample data is missing.
Crossvalidation allows one to find regions where data is over/underestimated and to compare the validity of different geostatistical approaches for data interpolation. By comparing sample values and their estimations, one can reveal data points having the greatest standardized errors. Such values demand additional attention because they may be data outliers. The result of validation and crossvalidation techniques can be analyzed further by using graphics tools and spatial analysis functions.
Some additional tools are available for the geostatistical data analysis. In this article we used declustering techniques, which allows one to improve estimation if measurements are sampled preferentially (see for example Deutsch, Journel, 1998).
In the article we used the following geostatistical approaches: simple, ordinary, indicator, probability and disjunctive krigings as well as ordinary and indicator kriging modifications for binomial data.
Simple (in this method the estimation of the mean can be established a priori based upon a different data set from the data used for the present estimation; for example, declustering mean) and ordinary (with unknown mean) krigings (Gandin, 1963) can be applied successfully to data having, or close to, a multivariate Gaussian distribution. The assumption of the multivariate Gaussian distribution is rarely realized in practice. To check for the multivariate normal distribution of the data, one should first transform the initial data to a univariate normal distribution (i.e. performing linearization of the original data). Such transformation performed prior to the variogram/covariance analysis allows one to reduce the variability in the original data, making the variogram modeling more reliable and stable. The normal score transform function can be defined for any continuous cumulative distribution function. Transformation of the distribution of the initial data set into univariate Gaussian distribution does not guarantee that all of the initial data will be transformed into an exact multivariate Gaussian distribution. It is, however, relatively easy to check for bivariate normality of the transformation (Deutsch, Journel, 1998). If bivariate distribution is fulfilled, transformation could be adopted. Simple and ordinary kriging provide an optimal estimation in the class of linear models. Suggestion about normality of the transformed data allows one to calculate the confidence interval for the estimation and to present results of the estimation as a map with the desirable level of significance and as a map of the probability that selected critical level is exceeded or not, in addition to the map of estimation.
Indicator kriging (Journel, 1988) is an example of the techniques in which the uncertainty model depends only on the information available. In indicator kriging a 01 transformation is introduced, which will be the new variable for which we will compute variograms and carry out kriging methods. Among the shortcomings of indicator kriging is the loss of information after coding data through indicator functions. For example, if the data values are in the range of 1 to 100 and the indicator value (cutoff) selected is 20, then the data points with values 21 and 99 will be interpreted as being equivalent. It is possible, therefore, in such situations to use a soft indicator function. This indicator function can be prepared by the user based upon knowledge of the data sets and processing by ordinary kriging. One of the possible solutions is to remove data near the threshold from consideration.
Probability (Journel, 1988) kriging is considered an improvement over the indicator kriging method in the sense that the data is used more completely for the task of estimating conditional probabilities that given thresholds are exceeded. This nonparametric approach uses the original data set in addition to indicators.
Disjunctive kriging (Rivoirard, 1994) is also based on a specific nonlinear transformation of the original data. The first step in this type of estimator is the normal score transformation of the data. The next step is to express the normalized function as an expansion of hermitian polynomials. It should be noted that a sufficient number of hermitian polynomials depends on searching neighborhoods and to receive reliable results of estimation we should use a different number of hermitian polynomials for different locations. It is also assumed that the condition of bivariate standard normal distribution is fulfilled. It should be noted that bivariate normal distribution is weaker than the multiGaussian condition. As mentioned above, it is possible to check the bivariate normality of the transformation. If bivariate normality exists, disjunctive kriging can produce better results than other geostatistical methods. If it is not so, it is desirable to use another type of kriging. Validation and crossvalidation techniques can be used to confirm this.
Environmental applications often include a collection of data, that can be processed by cokriging. Cokriging is the logical extension of kriging to situations where two or more variables are spatially interdependent and where the one whose values are to be estimated is not sampled as intensively as the others with which it is correlated. In other words, cokriging is used when one variable is of principal interest and the other variables are used to enhance the estimation of the first (i.e. the supplementary variable is used to compensate for a lack of data for the primary variable). Analysis of simple and crosssemivariograms, which are created during this process, allows one to check the spatial correlation between different variables. In general, using cokriging techniques will always result in more precise estimations. The influence of the secondary data set over the primary data set is higher if the primary data set includes the greatest error and if correlation between the data sets occurs. If the correlation between the data sets is negligible, we will receive the same result as kriging one variable. The difference of data estimation based on kriging and cokriging approaches could also help to find outliers.
LongLived Radionuclides Soil Contamination
Since the accident at Chernobyl, public safety in the territories with longterm radioactive contamination has assumed a fundamental importance. Investigation of soil radioactivity is a necessary step in the detailed description of the sources of radioactivity exposure because contaminated soil is the main supplier of radioisotopes into other components of the biosphere. This is the only possible way at the present time of determining the hazard from a radiators (mainly ^{238240}Pu and ^{241}Am) to human health since the direct dose measurements in this case are infeasible.
According to International Basic Standards for Protection against Ionizing Radiation and for the Safety of Radiation Sources (International…, 1996), radiation dose limits in unrestricted (uncontrolled) areas shall be such that an individual will not receive a dose to the whole body in excess of 1 millisievert (mSv) per year. The dose is the cumulative exposure to radiation actually received by an individual. In practice, however, classifications of the contaminated territories are based on the level of soil contamination and usually on ^{137}Cs soil contamination. (It is usually suggested that there exists a linear correlation between soil contamination and dose loads on population.) Comparison between disease rates and radiation exposure is usually also based on ^{137}Cs soil contamination. However, for Belarusian territories with less than 370 kBq/m^{2} (relatively low contaminated territory) there is no correlation between food and soil radiocesium contamination (i.e. the linear correlation between internal dose and soil contamination does not exist for most Belarussians [Konshin, 1992; Krivoruchko & Makeichik, 1997; Krivoruchko, 1997b]).
Figure 1a presents ordinary kriging estimation of ^{137}Cs
soil contamination for the Belarus territory based on 15,505 measurements
in 1992 (Krivoruchko, 1997a; Krivoruchko, Gribov, 1997b). Results of the
interpolation are presented as a hillshaded map and a threedimensional
surface. The arrow indicates the same location for these two surfaces.
a) b)
Figure 1. a) ^{137}Cs soil contamination for the Belarus territory, Ci/km^{2}. b) ^{137}Cs soil contamination near Volozhin City, Bq/m^{2}. Ordinary kriging estimation.
Figures 1b, 2a, and 2b present the ordinary kriging estimation
of ^{137}Cs soil contamination for three different clusters: the
western area between Minsk and Grodno, Figure 1b; the eastern area between
Gomel and Mogilev, Figure 2a; and the area between Gomel and Brest on the
south of the Republic (this is the area with the highest rates of radiationrelated
disease in Belarus), Figure 2b. Figure 2c presents the error of ordinary
kriging estimation for this area. In this figure, the symbol x represents
data collection locations (settlements). The darker the color in Figure
2c, the higher the error of estimations. Generally, errors depend on the
density of samples and on data variability in the locality. The highest
error of estimations is in the southern uninhabited part of this region,
which was a former Soviet army testing area.

b) c) 
Figure 2. ^{137}Cs soil contamination for different parts of Belarus, Bq/m^{2}. a) Eastern area. b) South central area. c) Error map for the south central area. Ordinary kriging estimation.
To prepare the maps for the different subregions of Belarus, we used
between 500 and 800 samples for each subregion and made kriging consider
anisotropy. Figure 3a presents variograms for the western area of radioactivity
in Belarus (set of the variogram in various directions) and Figure 3b the
eastern area (two directional variograms in perpendicular directions).
See the associated maps of contamination in Figures 1b and 2a.
It is clear that the structure of the contamination depends on the direction and distance from the Chernobyl and it is not appropriate to use standard kriging procedure for the entire data set. To overcome this problem it is possible to use either a socalled stratified kriging (i.e., doing interpolation for different subregions separately and then combining results of the interpolation on a single map) or an approach based on the combination of artificial neural networks and geostatistics to prepare and use a continuous field of correlation between data instead of using unique spatial correlation for the whole area.
The next example of radionuclide soil contamination shows the distribution of ^{90}Sr in 1994 (kBq/m^{2}) in southeast Belarus, in Gomel province (Krivoruchko, 1997a). The main parameters of the data are as follows: number of samples is 1486, minimum is 0, maximum is 97.3, median is 6.3, mean is 10.7, standard deviation is 12.6, skewness is 2.46, and kurtosis is 10.8.
This data set has an interesting feature. Measurements were not made in the most contaminated areas near Chernobyl and to the east of the Gomel province (see Figures 1  3 with data for the whole Belarus territory). People from these areas were evicted. It is important to see how different methods will interpolate the data in such areas.
The map in Figure 4a indicates the locations of sample
points by circles 5 km in diameter. We used different deterministic and
statistical methods for processing this data set, and parts of the maps
were published in (Krivoruchko, 1997a). As was expected, the best results
were obtained using different geostatistical approaches. Figure 4b illustrates
the successful interpolation of data in the most contaminated areas in
spite of missing sample points using probability kriging.
a) b)
Figure 4. ^{90}Sr soil contamination in the Gomel province of Belarus. a) Symbol map. b) Conditional probability that a threshold of 15 kBq per sq.m was exceeded using the probability kriging.
The geostatistical approach produces reliable results and allows one to know errors of estimations and, as a key advantage, one can produce maps of errors of estimations and maps of probability for a given variable to exceed a chosen threshold. The next example of using this extremely useful method for decision making is illustrated below using plutonium and americium soil contamination data.
Until the middle of the twentieth century plutonium was practically absent in the human. Its concentration in uranium and thorium ores is estimated to be about 10^{9}10^{7} Bq/g. Substantial amounts of plutonium were released into the natural environment after the beginning of ground nuclear weapons testing in the middle 1960s. In this period, plutonium concentrations in the upper soil layers of Belarus had reached the values of 10^{5}10^{4} Bq/g. After the accident at the Chernobyl NPP the concentrations of plutonium in soils in the southern areas of Belarus increased to 0.10.2 Bq/g. Unfortunately, the biogeochemical activity of plutonium remains poorly investigated.
A peculiarity of the plutonium deposition process after the Chernobyl accident was that the essential part of the radionuclide was associated with fuel particles (hot particles) of different size. This fact can partly explain the existence of outliers and a weaker correlation among the data in some subregions.
We have investigated the spatial structure of a
radiators in (Krivoruchko et al., 1998a; Krivoruchko et al., 1998b) and
present some of the results from these articles below. Figure 5a presents
the conditional probability that a threshold 370 Bq/m^{2} of ^{239,240}Pu
soil contamination in the Gomel province in 1992 was exceeded, calculated
using indicator kriging. Figure 5b presents the conditional probability
that a threshold 100 Bq/m^{2} of ^{241}Am soil contamination
was exceeded using probability kriging estimation.
a) b)
Figure 5. a) ^{239,240}Pu soil contamination in the Gomel province in 1992. Conditional probability that a threshold 370 Bq/m^{2} was exceeded using indicator kriging. b) ^{241}Am soil contamination. Conditional probability that a threshold 100 Bq/m^{2} was exceeded in 1992 using probability kriging.
Because maps of the distribution area of a radiators have not been published, the common position in the numerous articles devoted to studying the consequences of Chernobyl is to say that dangerous levels of contamination from plutonium are only in the restricted zone where people do not live and that the situation is under control. However, map in Figure 5b indicates that a large number of people are living in the danger zone now: cities, that are represented in Figure 5b are district centers in which population is more than 20,000. Gomel has about 300,000 inhabitants.
The radiation conditions within the restricted 30kilometer zone around Chernobyl have stabilized. Landscaperelated redistribution of radionuclides is completed and pine needles are completely renewed. The population left these territories during the first weeks after the accident. Since that time the restricted zone has been under thorough investigation on various aspects of the problem of radioactive territory contamination. Collecting and processing data on the concentration of longlived radionuclides in soils is of great importance for deciding on the acceptability of renewing agricultural activity and returning the formerly relocated population. The total amount of a radiators is changing over time due to the radioactive decay of relatively shortlived radioisotopes of ^{238}Pu to ^{234}U and ^{241}Pu to ^{241}Am, as well as due to vertical and horizontal migration of isotopes in soil. The migration velocity of a radiators in the most abundant podzol soils is very small. Over the next hundred years, these radionuclides will remain within the cultivated layer of soil, so we did not take the migration into account.
For the purpose of population dose assessment, it is necessary
to account for all a radiators, having many
common features of transfer and retention in tissues. The map in Figure
6a illustrates the predicted distribution of the total amount of a
radiators over Belarus in 2059, at the time of the maximum total concentration
of a radiators.
a) b)
Figure 6. a) Prediction of the a radiators in 2059. b) Areas of soil contamination with a radiators exceeding 740 Bq/sq.m in 1992 (solid pattern) and prediction for 2059 (dashed pattern). Ordinary kriging estimations.
Figure 6b shows territories with a radiators surface density of 740 Bq/sq and higher. Solid pattern reflects the contamination in 1992 and the dashed pattern is the prediction for 2059.
The sample area and the total amount of the activity of each radioisotope, as well as the sum of activities of the investigated a radiators in 1992 and 2059, were estimated using the Monte Carlo method. The sample area was estimated to be 1,180 km^{2}. The total activity of a radiators in the restricted Chernobyl zone in 2059 will exhibit almost a twofold increase of the activity compared to 1992 and the danger zone will expand. We believe that it will be impossible to live in the restricted 30kilometer zone around Chernobyl and its immediate borders for the rest of this century and into the next as well. More information about data, variography, correlation between a radiators and cokriging interpolation can be found in (Krivoruchko et al., 1998a; Krivoruchko et al., 1998b).
Spatial Distribution of RadiationInduced Thyroid Cancer
In this section of the article, spatial distribution of radiationinduced thyroid cancer in Belarus is presented based on the results of (Krivoruchko, 1997a; Krivoruchko, Naumov, 1997; Lutsko, Krivoruchko, 1996; Lutsko et al., 1997).
About five years after the Chernobyl nuclear accident, the first suspicions arose about an increased incidence of thyroid cancer in children (Kazakov et al., 1992; Prisyazhiuk et al., 1991). A debate started on whether this increase was real and could be attributed to the release of radiation as a consequence of the accident or whether it was an artifact. In the meantime, with the growing number of cases, the critical voices became silent. The cumulative incidence rate over ten years (1986  1995) reached as many as 17.2 cases per 10,000 children in the Bragin district to the east of Chernobyl. It is followed by the Narovlya (16.8 cases per 10,000 children) and the Hoiniki district (12.8 cases per 10,000 children). Last year in Belarus there were about 100 new cases of thyroid cancer in children. Prior to Chernobyl the rate was about one case per the entire Belarus child population.
The data on the distribution of thyroid cancer among the population of Belarus was provided by the Minsk Oncology Institute. The database includes 3,773 case records, 902 of which correspond to the preChernobyl period, and the remainder to postChernobyl. To define the incidence rate of thyroid cancer in districts and settlements we used information from the 1979 and 1989 censuses of population and information about district population from the annual demographic reports. The age structure of thyroid cancer morbidity in Belarus was investigated in numerous articles (Lutsko et al., 1996).
The cumulative incidence figures for each district were assigned to the coordinates of the district centroids. The 117 incidence values were then considered a spatially distributed variable, that could be interpolated with regard to the spatial autocorrelation found in the variogram.
The feature of the spatial problem of thyroid cancer spatial distribution is that the population strongly differs from one settlement (district) to another. The average number of thyroid cancer cases is two to three per 10,000 over the whole postChernobyl period. Thus, the rate of incidence could be correctly defined by the average value at the settlement with population more than 100,000. When population of a settlement is close to 100,000 or less, the rate of incidence becomes inaccurate (there are about 30 districts having less than 30,000 inhabitants). Some cases relate to villages having twenty or more people. When considering the problem of morbidity per settlement one must take into account that uncertainty will rise. In such cases the direct use of the average rate of incidence would be incorrect. One must take into account the initial data structure, which is binomial.
For spatial interpolation of thyroid cancer disease we used a modification of ordinary kriging, which was originally proposed (McNeil, 1991). In order to create a map of the conditional probability that a predefined level is exceeded, we proposed soft indicator kriging modification for binomial data. It excludes from consideration some of the values lying in close vicinity to the threshold value (Krivoruchko, Gribov, 1997a). Thus, in mapping the risk of thyroid cancer development we considered both the distances between settlements or centers of districts and the number of inhabitants. We can say that large settlements with large populations will carry the most accurate characteristics of the epidemic situation in their vicinity, but nevertheless information related to small settlements will not be lost.
Figure 7a represents the map of the estimated risk of
thyroid cancer in the postChernobyl period (2,756 cases of disease) among
all of the population. Figure 7b presents a map of the estimation of thyroid
cancer risk for children from 1986 to 1995. The analysis was limited to
528 individuals who were under 15 in 1986. It is presented as the specific
probabilities that a predefined threshold of one case of the disease per
10,000 children was exceeded.
a) b)
Figure 7. a) The estimated risk of thyroid cancer in the postChernobyl period (2,756 cases). Ordinary kriging for binomial data. b) Conditional probabilities that the risk of thyroid cancer in children is higher that one per 10,000. Soft indicator kriging estimation.
Prior to Chernobyl, thyroid cancer incidence was greater in the northern part of Belarus, attaining the highest levels in the Vitebsk district. After the Chernobyl accident the pattern of thyroid cancer incidence substantially changed. Besides the area near Vitebsk, many new areas showed an increased incidence of disease. The territories located along the western trace of the radioactive release suffered the greatest. New areas of increased thyroid cancer risk emerged around Minsk, though iodine doses in this area were relatively low.
There is strong evidence that the increased incidence of childhood thyroid cancer is due to radiation exposure as a result of the Chernobyl accident, based on the geographical and temporal distribution of the cases. The map of the estimated risk of thyroid cancer agrees with the map of the reconstructed surface iodine pollution of the Belarus territory (Krivoruchko, Naumov, 1997; Lutsko et al., 1997) as well as with dose loads on the population in the first days after the accident (see next section).
Dose Loads on Population in the First Days After the Chernobyl Accident
Due to the unfortunate fact that the measuring equipment was inadequate to properly monitor the scale of radiation exposure during the early period of the accident, detailed direct information on the deposition of the shortlived radionuclides and the doses to the population has been irretrievably lost. Now the only way to reconstruct the dynamics of the radioecologic situation of the initial period of the Chernobyl accident is to make a retrospective assessment of radiation exposures related to the shortlived radionuclides. This can be done by processing extensive empirical information on radiology, meteorology, and epidemiology, and through comparison and validation of the data.
Radioisotopes from the nuclear plant released from April 26 to May 5, 1986, were transferred and deposited over Belarus in a substantially irregular pattern. For the purpose of dose assessment, one must consider the daily structure of radioactive releases in different directions, starting at April 26. For example, over the territories to the west of Pripyat, the former residential community of the Chernobyl nuclear plant, the main clusters of the shortlived radionuclide contamination were formed during the period of April 27  30, 1986. To correctly assess the contribution made by shortlived radionuclides to doses received by the population, it is necessary to consider numerous processes, starting with the spectrum of the released activity. To prepare the data for further analysis using GIS and geostatistics, we used 114 exposure dose measurements made on May 10, 1986, and took into account the following (Krivoruchko, Naumov, 1997):
Figure 8 presents the results of the reconstruction of the thyroid dose loads from ^{131}I, based on consumption of contaminated milk by the adult population between April 26, 1986, and May 5, 1986. It was suggested that there was no migration at this time and that people used local foods. Tragically, it is a very reliable suggestion since evacuation from Belarus commenced on May 1, 1986, and by May 5, 1986 only 11,035 people had been relocated.
For the interpolation we used simple cokriging estimation
with the normal score transformation of the primary data. The secondary
data set was cumulative data of thyroid cancer incidence in children from
1986 to 1995 (117 data samples, one for each of the districts in Belarus).
a) b)
Figure 8. Thyroid dose loads (in mSv) in adult population from April 26, 1986, to May 5, 1986. Simple cokriging estimation. Secondary data set is thyroid cancer incidence in children from 1986 to 1995. a) External gamma doses. b) Internal doses.
We used ordinary kriging for estimating and mapping the external and internal doses of the population (Krivoruchko, Naumov, 1997). The cokriging approach allows one to reduce the error of estimation when the variables used are related by a common process, as in our case.
It should be noted that the distribution of shortlived radionuclides is different from longlived radionuclides (see the section "LongLived Radionuclides Soil Contamination") as well as the patterns of dose on population in the very first days of the Chernobyl accident and for 19871998 (Krivoruchko, 1998).
^{137}Cs Food Contamination in 1993
Currently, the main sources of radiation hazard to the population living in the contaminated regions are internal exposure from food and external exposure from gamma dose rate in the air. This section is based on the measurements of ^{137}Cs contents in 83 types of food, that were carried out in the Byelorussian Institute of Radiation Safety in 1993. The cases of ^{137}Cs content in food exceeding the upper permissible level were published in information bulletins (Nesterenko, 1996), which allows one to identify families under high radiation risk. The initial database containing 53,207 records on ^{137}Cs concentration in 83 types of food was available for the present investigation. Radiocaesium contamination is distributed very nonuniformly both geographically and within different types of food (see Table I).
Table I. Exceeding the ^{137}Cs upper permissible
level in some types of food in rural settlements of Mogilev, Gomel, and
Brest provinces in 1993.
Food  UPL, Bq/kg  Number of Settlements  Number of Measurements  Exceeding UPL, % 
Bilberry  185  317  1,383  61.03 
Mushrooms  370  292  1,123  56.0 
Milk  111  675  19,111  14.85 
Pork  185  229  234  14.10 
Sour cream  111  83  242  12.81 
Well water  18.5  243  2141  8.78 
Carrots  185  252  1,439  5.84 
Cabbage  185  182  590  4.41 
Potatoes  370  472  4,996  1.64 
For decision making it is important to find regions where
it is dangerous to consume food. We have prepared probability maps for
each type of food based on the probability kriging approach. Figure 9 presents
maps of conditional probabilities that the upper permissible levels for
mushrooms and bilberries were exceeded.
A number of factors influence the uptake of radionuclides from soil to plants including the level of soil contamination, the soil type, and the type and extent of countermeasures. We have analyzed the spatial correlation between ^{137}Cs in the food and in the soil and found that this dependence is very complicated.
^{137}Cs milk contamination was changed significantly
from 1988 to 1993. Figure 10 presents results of conditional probability
estimation that levels of 37 Bq/kg in 1993 and 185 Bq/kg in 1988 were exceeded.
a) b)
Figure 10. ^{137}Cs concentration in milk. a) Conditional probability that level 37 Bq/kg was exceeded in 1993. b) Conditional probability that level 185 Bq/kg was exceeded in 1988.
The maps presented in Figures 9 and 10a are different from the map of ^{137}Cs soil contamination (Figures 1 and 2), especially in relation to the western passage of the radioactive cloud in the very first days after the Chernobyl accident.
The estimation of the internal exposure was made based on the following food components: milk, potatoes, vegetables, fruits and berries, meat, mushrooms, fish, and bread. We used the dietary habits of the adult rural population of Belarus and information about the diet of children from the Chernobyl zone (Nesterenko, 1996).
Figure 11a presents the structure of the internal dose
loads component from food intake for the adult population in 1993, that
was based on the information for 120 rural settlements for which we had
at least 50 measurements for each food component.
The high radiation risk area covers almost all of the southern part of Belarus. Visual analysis of maps allows one to conclude that the most contaminated food is produced in the rural settlements around the Stolin district, which differs from the official point of view. As was mentioned in (Nesterenko, 1996), practically all children in this area have health related problems. More information about ^{137}Cs food contamination can be found in (Krivoruchko, 1997b; Krivoruchko, Makeichik, 1997; Krivoruchko, 1998).
Conclusion
Many methods of spatial data interpolation exist that can analyze data in a reasonable time, however, the choice of a particular method is critical if data is sparse or includes errors. In the section Spatial Distribution of RadiationInduced Thyroid Cancer, we presented the results of estimating thyroid cancer incidence for the adult population of Belarus after the Chernobyl accident (see Figure 7a). Let us consider the environmental data processing using this example based on some other commonly used methods of spatial data interpolation.
Interpolated values at unvisited points in Figure 12a were obtained by one of the most commonly (due to its simplicity) used method  inverse squared distance weighting. Inverse squared distance weighting interpolation is defined by the following formula , where parameter p equals 2. The surfaces resulting from a weighted average interpolation depend on the parameter p and on the size of a window from which the sample data was drawn. In the case when p» 20, the map will be very close to a Voronoi map. The inverse squared distance weighting approach to data interpolation sends the user in search of the "best" parameters, which, however, cannot be based on sound mathematical considerations.
Figure 12b presents the results of a trend surface analysis
technique global polynomial interpolation. This simple deterministic method
allows one to produce smooth maps that depend on the order of the polynomial.
a) b)
Figure 12. Risk of thyroid cancer in the postChernobyl period among all Belarus population. a) Inverse square distance weighting method. b) Global polynomial interpolation. Twodimensional surface of 5order.
Trend surface interpolation is highly sensitive to outliers (extremely high and low values) and it is not an exact interpolator (i.e. resulting isolines do not pass through each of the measurement points).
It is possible to improve this method by using data in
the specified neighborhood. Additional improvement of this method can be
achieved using the weighting coefficients of the polynomial (Gandin, 1963).
Figure 13a displays the results of the local polynomial interpolation by
a thirdorder polynomial. Figure 13b presents the results of standard indicator
kriging estimation of the conditional probability that the risk of thyroid
cancer is higher than 25 per 100,000. The difference between the map in
Figure 13b and the maps in Figures 12 and 13a is apparent.
a) b)
Figure 13. Risk of thyroid cancer in the postChernobyl period among all Belarus population. a) Thirdorder local polynomial interpolation. b) Conditional probability that risk of thyroid cancer is higher that 25 per 100,000. Indicator kriging estimation.
The construction of the maps in Figures 12 and 13a raise some important questions, that include:
Acknowledgments
I thank Hugh Keegan and Peter Petri for useful discussions.
References
Deutsch C.V., Journel A.G. (1998) GSLIB: Geostatistical Software Library and User’s Guide. Oxford University Press, New York.
Gandin L.S. (1963) Translation: Objective Analysis of Meteorological Fields. Israel Program for Scientific Translations. Jerusalem, 1965.
International Basic Standards for Protection against Ionizing Radiation and for the Safety of Radiation Sources. International Atomic Energy Agency. Vienna (1996).
Journel A.G. (1988) Nonparametric Geostatistics for Risk and Additional Sampling Assessment. In Lawrence H.Keith (ed), Principles of Environmental Sampling. American Chemical Society, pp. 4572.
Kazakov V.S., Demidchik E.P., Astakhova L.N. (1992) Thyroid Cancer after Chernobyl. Nature, 359:21.
Konshin O.V. (1992) Transfer of Cs137 From Soil to Grass: Analysis of Possible Sources of Uncertainty. Health Physics, V63, pp. 307315.
Krivoruchko K. (1997a) Geostatistical Picturing of Chernobyl Fallout and Estimation of Cancer Risk Among Belarus Population. Third Joint European Conference and Exhibition on Geographical Information, Vienna, Austria, April 1997, Vol. 2, pp. 676685.
Krivoruchko K. (1997b) Spatial Structure of Food Contamination with ^{137}Cs and Estimation of Longterm Internal Dose Loads on Population of Belarus. International Conference on Low Doses of Ionizing Radiation: Biological Effects and Regulatory Control. Seville, Spain, 1721 November 1997, pp. 152155.
Krivoruchko K. (1998) Radionuclides Soil Contamination, Dose Estimation and Dynamic of Cancer Incidence in Belarus. Minsk, in print (in Russian).
Krivoruchko K., Gribov A. (1997a) Estimating the Risk of Cancer in Belarus: Lung Cancer Example. Ecological Antropology, vol. 2, MinskLublinLodz, pp. 423427.
Krivoruchko K., Gribov A. (1997b) The Analysis of Spatial Structure and Variability of the Chernobyl Fallouts by Geostatistics. International Symposium "Actual Problems of Dosimetry", Minsk, Belarus, 2830 October 1997, pp. 183194 (in Russian).
Krivoruchko K., Makeichik A. (1997) Spatial analysis of radiocaesium food contamination in rural settlements of Belarus. International Symposium "Actual Problems of Dosimetry", Minsk, Belarus, 2830 October 1997, pp. 195199.
Krivoruchko K., Gribov A., Figurin I., Karebo S., Pavlushko D., Remesov D., Zhigimont A. (1997) Mapstudio: The Specialized GIS Integrating Possibilities of Geostatistics. Third Joint European Conference and Exhibition on Geographical Information, Vienna, Austria, April 1997, pp. 187196.
Krivoruchko K., Mironov V., Krivoruchko T., Kudrjashov V. (1998a) About the Plutonium Contamination of the Belarus Territory. Ecological Anthropology, Proceedings of the International Conference: Human Ecology in the PostChernobyl Period—Children of Chernobyl, Minsk, pp. 320323.
Krivoruchko K., Mironov V., Gribov A., Krivoruchko T., Korolev V., Kudrjashov V. ^{238}Pu, ^{239,240}Pu and ^{241}Am Soil Contamination in the Restricted Chernobyl Zone. (1998b) Ecological Anthropology, Proceedings of the International Conference: Human Ecology in the PostChernobyl Period—Children of Chernobyl, Minsk, pp. 323325.
Krivoruchko K., Naumov A. (1997) Reconstruction of dose loads on population in the initial period of the Chernobyl accident and estimation of thyroid cancer risk in Belarus. International Symposium "Actual Problems of Dosimetry", Minsk, Belarus, 2830 October 1997, pp. 167173.
Lutsko A., Krivoruchko K. (1996) Analysis of methods and reconstruction of iodine doses from the Chernobyl release in Belarus. International Conference "One Decade After Chernobyl". Vienna, Austria, 812 April 1996, Vol. 2.
Lutsko A., Milutin A., Krivoruchko K. (1996) Health Characteristics of Children and Adults in Belarus, which Suffered from the Chernobyl Accident. 10 Anni da Chenobyl: Ricerche in Radiologia, Monitoraggio Ambientale e Radioprotezione. Trieste, marzo 1996, pp. 2938.
Lutsko A., Krivoruchko K., Gribov A. (1997) Analysis of the Methods and Reconstruction of Iodine Doses from the Chernobyl Release in Belarus. International Symposium "Actual Problems of Dosimetry", Minsk, Belarus, 2830 October 1997, pp. 512.
Maignan M., Krivoruchko K. (1997) Integration GIS and Geostatistics: A Software and a Case Study. Proceedings 18^{th} ICA/ACI International Cartographic Conference. Stockholm, Sweden, 2327 June 1997, Vol.2, pp. 925932.
McNeil L. (1991) Interpolation and smoothing of binomial data for the southern african bird atlas project. South African Statistical Journal. V.25, pp. 129136.
Nesterenko V.B. (1996) Radioecological Monitoring of Food in Rural Settlements in Belarus. Pravo i economika. N 19, Minsk, 19961998. In Russian.
Prisyazhiuk A, Pjatak O, Buzanov V.A, Reeves G.K, Beral V. (1991) Cancer in Ukraine, postChernobyl. Lancet, 338:13345.
Rivoirard J. (1994) Introduction to Disjunctive Kriging and NonLinear Geostatistics. Oxford Univ. Press. New York.
Sakharov A. (1958) Radioactive Carbon from Nuclear Explosions and Nontreshold Biological Effects. Soviet Journal of Atomic Energy, vol.4, N6.
In May 1998, I began a three year contract with ESRI for the Software Development Team. Prior to that, I was an Assistant Professor and Director of the GIS Laboratory at the Sakharov Institute of Radioecology in Minsk, Belarus. My responsibilities include the development and utilization of the GIS Software MapStudio, encompassing databases from the Belarussian National Cancer Registry; data pertaining to radionuclides contamination in soil and food from Belarussian National Academy of Sciences; the development of GIS/spatial statistics curriculum; and supervision of Ph.D. and graduate candidate research projects pertaining to GIS applications and spatial data analysis based on geostatistical approach. Prior to joining the Sakharov Institute, I was a divisional chief at the Belarussian National Academy of Sciences, Nuclear Power Institute, focusing on issues concerning radioecology and epidemiology.