
Olivier Perrin
INRA, Domaine SaintPaul, Site Agroparc, 84914 Avignon Cedex 9 France
Email: perrin@avignon.inra.fr
Serge Iovleff
Laboratories SABRES, IUP VannesTohannic, rue Yves Mainguy, 56000 Vannes France
Email: Serge.Iovleff@univubs.fr
In the last decade, a useful model for nonstationary random fields has been developed. This consists of reducing the random field of interest to isotropy via a bijective bicontinuous deformation of the index space. Then the problem consists of estimating this space deformation. We propose to estimate this space deformation using a constrained continuous version of the simulated annealing for a Metropolis dynamic.
Assumption of isotropy is clearly violated for many, if not most, spatial environmental phenomena. Factors such as topography, local pollutant emissions, and meteorological influences may cause such assumptions to be violated. This has led to research into modelling spatially nonstationary second order structure, as reviewed in [7]. For dealing with this nonstationarity, Sampson and Guttorp [19] have developed a model that consists of reducing the correlation function r(x,y) of the spatial phenomenon of interest, modelled by a random field , to an isotropic one as follows
(1) 
where denotes the Euclidean norm in , represents a bijective deformation of the geographic coordinate system and is an isotropic correlation function. In the sequel, we refer to the geographical coordinate system as the Gspace , where G stands for geographical, and to the deformed coordinate representation as the Dspace , where D stands for deformed.
Unlike the classical geostatistical models, nonstationarity through second order moments is thus taken into account and model (1) gives opportunity to enlarge the class of models for studying spatial environmental random fields.
When is strictly decreasing, Perrin and Meiring [16] prove the uniqueness of both the deformation and up to a homothetic Euclidean motion for and up to a scaling for . Perrin and Senoussi [18] give the general form of the deformation that reduces a nonstationary random field in the way (1), under smoothness assumptions.
Our concern in this paper is the estimation of the space deformation in model (1) when the isotropic correlation depends on a parameter , . From now on, stands for the isotropic correlation function. This estimation is based on T repetitions (independent and identically distributed observations) of Z at each of n distinct monitoring sites in G, which may be irregularly located; we suppose that G represents the convex hull of these sites. We denote these records using Z_{t}(x_{i}), , .
So far, the deformation has mostly been estimated with the help of a nonparametric approach. This nonparametric estimation of has already been extensively developed, [11], [12] and [13], and applied in analyses of solar radiation [19], acid precipitation [8], rainfall [14], air pollution [3], tropospheric ozone [20], ...This approach consists of modelling using a pair of thinplate splines and minimising a penalised weighted least squares criterion to estimate the Dspace coordinates of the monitoring sites x_{i}, . This treatment is carried out with a Marquardttype algorithm. At least, two drawbacks are associated with this method: (i) bijection condition is not ensured; (ii) the fitting of the model becomes a challenging numerical problem with dimensionality roughly proportional to the number of fixed monitoring sites [12].
Perrin and Monestiez [17] propose a parametric approach. They model using a composition of a "small'' number of bijective elementary radial basis deformations. Then they use a least squares criterion to estimate the parameters; this criterion is minimised in a classical way (Marquardttype algorithm). Let us precise two disadvantages of this method: (i) so far there is no rational choice of the relevant number of elementary deformations; (ii) there is a considerable dependence on the starting values in the minimisation procedure.
To avoid all these disadvantages, we propose here to estimate in a first step the Dspace coordinates , , by minimising an objective function with a stochastic algorithm, the continuous version of the simulated annealing for a Metropolis dynamic. To ensure a bijective correspondence between the n points in the Gspace and their estimated deformations in the Dspace, we impose some nonfolding constraints in the algorithm. In a second step, to estimate the deformation in the whole Gspace, we use a piecewise affine interpolation of the points x_{i} in the Gspace and the estimations of their deformations in the Dspace, .
The paper is structured as follows. In Section 2, we address the estimation problem of the deformation with the simulated annealing algorithm together with some nonfolding constraints. In Section 3, we illustrate our nonparametric approach with precipitation data from 20 sites in the LanguedocRoussillon region of France. Section 4 points out one possible application where the nonstationary model of a correlation function using spacedeformation may be useful. More precisely, we give one idea of how spatial prediction should proceed in the new coordinate space and propose a cross validation study to demonstrate the possible improvement due to the deformation in predictions. Further, we apply this cross validation study to our precipitation data set. Finally, Section 5 discusses two extensions and future directions related to this work.
The repetitions of Z allows us to define the sample correlation estimates for each couple of sites (x_{i},x_{j}),
(2) 
where are the sample mean estimates: We denote using , , the positions of the sites in the Dspace and we define and as the estimations of and which minimise the following objective function
(3) 
with respect to the parameters and , and subject to some nonfolding constraints described in Paragraph 2.2.
We know from Perrin and Senoussi [18] and Perrin and Meiring [16] that the bijective deformation and the isotropic correlation in (1) are unique up to a homothetic Euclidean motion for the deformation and up to a scaling for the correlation. The practical implication of this result is that we can fix two site positions in the minimisation procedure to fix translations, rotations and scalings.
The objective function (3) is minimised using a continuous state version of the simulated annealing algorithm ([6], [10], [1], [5] and [2] for seminal references) subject to some nonfolding constraints.
Among others, simulated annealing has the following advantages:
For all these reasons and because of the high complexity of our objective function (3), we use simulated annealing for minimising (3). However, let us precise that in our particular context, the main advantage of this algorithm consists of its great flexibility. By this we mean that it can include nonfolding constraints to ensure a bijective estimation of .
We describe hereafter the different steps to minimise (3) with nonfolding constraints
 step 0: we start from the geographical configuration of the sites (Gspace). We set . Then to give an initial value to , minimise the objective function (3) with respect to (y(0) is hold fixed) with a Marquardtype algorithm. Let be the estimation of at step 0. Finally, take a sequence of "temperatures'' decreasing to 0 by steps of length n
 step k: let be the estimated position of the sites. The change proposition is: fix , choose only one site (candidate point) j uniformly among the n sites and move it locally at a position y with natural nonfolding constraints we precise hereafter. The other sites are hold fixed. Then:
These constraints account for global and local nonfoldings.
First, we embed G in a rectangle . We precise the function of this rectangle further.
Second, we construct the Delaunay triangulation (cf. [15] for instance) for the n geographical sites plus the 4 vertices of as well as the 4 midpoints of its 4 edges. For computing the Delaunay triangulation we use a program written by Shewchuk [21].
For each site x_{i}, , we identify all the triangles for which this site is a vertex. Then consider the polygon composed with the aggregation of these triangles and denote using its vertices, where q_{i} is the number of triangles. We assume that these vertices are ordered clockwise and denote using (A^{1}_{i}_{,l},A^{2}_{i}_{,l}) their coordinates , .
At step k of our algorithm, a point y_{j}(k1) is moved at a position y. To ensure nonfolding constraints, in other words to avoid that triangles overlap, we impose that this new position is chosen uniformly in a convex polygon which is the set of points M_{j}=(M_{j}^{1},M_{j}^{2}) such that:
with the convention q_{j}+1=1. This inequality means that for any M_{j} in the segment [y_{j}(k1),M_{j}] do not intersect any of the line (A_{j}_{,l},A_{j}_{,l+1}). The region corresponds to the marked area of Figure 1
Figure 1: The marked area corresponds to the nonfolding constraints.

These constraints mean that we impose moves in the simulated annealing algorithm that preserve the topological structure of the Delaunay triangulation the same before and after the deformation.
The rectangle allows us to write the previous constraints in a similar way for all the sites. Indeed, each site of G is included in a polygon for which the vertices are also vertices of the Delaunay triangulation for the n sites and the 8 points of . Without this embedding, we would have to distinguish between the boundary points of G (the vertices of the convex hull) and its interior points. Moreover, the embedding of G in prevents us from global folding of G in itself. In conclusion, we can say that the Delaunay triangulation accounts for local nonfoldings and that the embedding of G in accounts for global nonfoldings.
Any point x in G belongs to a (unique) triangle for which the vertices are three monitoring sites, say x_{i}, x_{j} and x_{k}. Thus, any point x is uniquely defined by its barycentric coordinates in the triangle as follows
x=a_{x}_{,i}x_{i}+a_{x}_{,j}x_{j}+a_{x}_{,k}x_{k}
where such that a_{x}_{,i}+a_{x}_{,j}+a_{x}_{,k}=1. We define the estimation at point x of the deformation as follows
where , and are estimated with the simulated annealing algorithm. So defined, the estimation holds the following features: continuous; bijective; an interpolant, i.e. , .
They consist of monitored precipitation data from n=20 sites in the LanguedocRoussillon, south of France, for which the geographical configuration is shown in Figure 2.
Figure 2: Site locations with department and coast contours at the top, and the corresponding Delaunay triangulation below.

These data are in the form of 10day aggregates, giving 6 records during November and December each year from 1975 through 1992 (T=108). These data hold the following features: (i) a low frequency of missing values; (ii) a similar altitude (between 0 and 200 meters) for all the sites so that an altitude correction is pointless; (iii) a homogeneous period in the year so that a seasonal adjustment is not necessary. As pointed out by Meiring [11], means and variances are positively related, and therefore, sample correlations between observations are calculated on the scale, after adding 1 to all observations.
The sample correlation estimates for a pair of sites is based on all the time points for which both sites has observations. The sample correlations for different pairs of sites are sometimes based on different numbers of observations, so that the sample correlations matrix is not positive definite; however, the model (1) fitted to the sample correlation remains positive definite.
The upper righthand plot in Figure 5 represents the geographical intersite distances versus the corresponding correlations 1<i<j<n defined by (2), to which one would fit an isotropic correlation model if we assume weak isotropy for Z (note this means that is the identity function in (1)). According to this plot, we decide to use the exponential model
so that the objective function (3) is rewritten as follows
(4) 
For minimising (4) we apply our constrained procedure described in Paragraph 2.2. The Delaunay triangulation for the 20 sites is shown in Figure 2. In the cooling schedule we take c_{0}=1000 and . Figure 3 shows the estimated deformation of the Delaunay triangulation.
Figure 3: Deformation of the Delaunay triangulation.

The estimated deformations of the Gspace coordinates also are illustrated by Figures 4. In this figure, the motion from one site in the Gspace to its estimated position in the Dspace is represented by an arrow.
Figure 4: Plot of the estimated spatial deformation from the Gspace (circles) to the Dspace (black dots).

The upper righthand plot of Figure 5 shows the fitted exponential correlation as a function of the Dspace coordinates. The Gspace has been stretched in regions of relatively lower correlations, and shrunk in regions of relatively higher correlations, so that the exponential correlation better models the correlations in the Dspace representation (minimum of the objective function is equal to 0.022) than in the Gspace representation (minimum of the objective function is equal to 0.155). Improvement of this fitting is illustrated by the interquartile intervals: the empirical correlations are less scattered in the Dspace than in the Gspace.
Figure 5: At the top and on the left, intersite distances in the Gspace versus empirical correlations and the fitted exponential correlation model (solid line). At the top and on the right, intersite distances in the estimated Dspace versus empirical correlations and the fitted exponential correlation model (solid line). At the bottom and on the left, boxplots of the empirical correlations as a function of the distances in the Gspace. At the bottom and on the right, boxplots of the empirical correlations as a function of the distances in the Dspace.

We point out one possible application where the space deformation model (1) may be useful. This is the prediction using Kriging. We refer to Cressie [4] for a presentation of the Kriging method.
In practice the mean and the variance fields are very seldom known, and must be predicted. As Høst et al., illustrate [9], separate modelling of the mean, variance and residual fields from monitoring data collected in space and time, may give very valuable information about the standard errors in spatial interpolation. Prediction of each of the mean, variance, and residual field contributes to the overall spatial interpolation errors; however, estimation of the mean and variance is not the topic of this work and we concentrate only on the prediction of the centred and standardised random field Z.
To predict the centred and standardised random process Z at any location we use the simple Kriging predictor
with where denote the vector of the Kriging coefficients, the correlation matrix of the vector and the correlation vector between the ``known sites'' and the ``unknown site'' Z(x).
One way to check whether the deformation may improve the prediction is a cross validation study. More precisely, at each time t, we set aside a site and we predict the value of Z at this site with the n1 remaining sites using the simple Kriging predictor described above. We repeat this operation for each of the n sites T times. Then we calculate the mean square error prediction that has to be compared with the one we would obtain using the fitted isotropic correlation function if we assume weak isotropy for Z. We apply this cross validation study to the precipitation data set for which we only keep the repetitions that contain no missing values (that is to say 92 repetitions). We have the following results
This means that in average the prediction at the monitoring sites is better in the Dspace than in the Gspace, so that we can claim that the model with deformation for the correlation of Z is better that the model without deformation. This has to be confirmed with the help of statistical tests to check if this improvement is significant or not.
Simulated annealing appears to be a suitable tool for estimating a nonstationary structure. Combined with nonfolding constraints, it gives the opportunity for the first time to estimate nonparametrically the bijective spatial deformation , in model (1), in such a way that the estimation remains bijective; however we wish to note two research questions
Olivier Perrin was supported by INRA, France; and by the European Union's research network "Statistical and Computational Methods for the Analysis of Spatial Data," Grant #ERBFMRXCT960095, while he was visiting the department of Mathematical Statistics at Chalmers University of Technology and Göteborg University in Sweden. The authors are grateful to Dominique Courault from Unité de Bioclimatologie, INRAAvignon, France, for providing the precipitation data set, and to Jonathan Richard Shewchuk [21] for providing the program that performs the Delaunay triangulations. We also are grateful to Xavier Guyon for giving us the idea of this work.
Arts, E. and T.J. Kors, 1989. Simulated annealing and Boltzman machines: stochastic approaches to combinatorial optimization and neural computing. Wiley.
Azencott, R., Ed., 1992. Simulated annealing: parallelization techniques. Wiley.
Brown, P.J., N.D. L.E. J.V. Zidek, 1994. Multivariate spatial interpolation and exposure to air pollutants, The Canadian Journal of Statistics 22, pp. 489505.
Cressie, N.A.C., 1993. Statistics for spatial data. Revised Edition: WileyInterscience, John Wiley and Sons, Inc. p. 900.
Geman, D., 1990. Random fields and inverse problem in imaging, L.N.M. No. 1427, Springer.
Geman, S. and D. Geman, 1984. Stochastic relaxation, Gibbs distribution, and the Bayesian restoration of images. IEEE Trans. Pattern Anal. Machine Intell. 6, pp. 721741.
Guttorp, P. and P.D. Sampson, 1994. Methods for estimating heterogeneous spatial covariance functions with environmental applications, in: Patil, G.P. and Rao, C.R. (eds), Handbook of Statistics XII: Environmental Statistics, Elsevier/North Holland, New York, pp. 663690.
Guttorp, P., P.D. Sampson, and K. Newman, 1992. Nonparametric estimation of nonstationary spatial covariance structure with application to monitoring network design, in: Walden, A. and Guttorp, P. Eds., Statistics in Environmental and Earth Sciences, Edward Arnold, London, pp. 3951.
Host, G., H. Omre and P. Switzer, 1995. Spatial interpolation errors for monitoring data, Journal of the American Statistical Association 90, pp. 853861.
Kirkpatrick, S., C.D. Gelatt, Jr. and M.P. Vecchi, 1983. Optimization by simulated annealing, Science 220, pp. 671679.
Meiring, W., 1995. Estimation of heterogeneous spacetime covariance, PhD thesis. University of Washington, Seattle, WA.
Meiring, W., P. Guttorp and P.D. Sampson, 1998.Computational issues in fitting spatial deformation models for heterogeneous spatial correlation, Computing Science and Statistics 29(1), (Scott, D.W., Ed.), pp. 409417.
Meiring, W., P. Monestiez P.D. Sampson and P. Guttorp, 1997. Developments in the modelling of nonstationary spatial covariance structure from spacetime monitoring data, in: Baafi, E.Y. and Schofield, N. (eds), Geostatistics Wollongong '96, Vol.1, Kluwer Academic Publishers, pp. 162173.
Monestiez, P., P.D. Sampson and P. Guttorp, 1993. Modelling of heterogeneous spatial correlation structure by spatial deformation, Cahiers de Geostatistique, Fascicule 3, Compte Rendu des Journées de Géostatistique, 2526 May 1993, Fontainebleau. Published by the École Nationale Supérieure des Mines de Paris.
Okabe, A., B. Boots and K. Sugihara, 1922. Spatial tessellations. Concepts and application of Voronoi diagrams. Wiley, Chichester.
Perrin, O. and W. Meiring, 1998. Identifiability for nonstationary spatial structure, Submitted for publication.
Perrin, O. and P. Monestiez, 1998. Modelling of nonstationary spatial structure using parametric radial basis deformations, to appear in: Soares, A., GómezHernandez, J. and Froidevaux, R. (eds), geoENV98, Kluwer Academic Publishers.
Perrin, O. and R. Senoussi, 1998. Reducing nonstationary random fields to stationarity or isotropy using a space deformation, Submitted for publication.
Sampson, P. and P.D. Guttorp, 1992. Nonparametric estimation of nonstationary spatial covariance structure. JASA 87, pp. 108119.
Sampson, P.D., P. Guttoro and W. Meiring, 1994. Spatiotemporal analysis of regional ozone data for operational evaluation of an air quality model, Proceedings of the Section on Statistics and the Environment, American Statistical Association.
Shewchuk, J. R., 1996. Triangle: engineering a 2D quality mesh generator and Delaunay triangulator, First Workshop on Applied Computational Geometry (Philadelphia, Pennsylvania), pp. 124133, ACM. Available online from http://www.cs.cmu.edu/ quake/triangle.research.html.