In geostatistics, the presence of outlying data is more the rule than the exception. Moreover, the statistical analysis of data contaminated by outliers requires caution, particularly when a spatial dependence exists. In order to take into account these possible outliers during the adjustment of the spatial process, a new modeling tool, called the substitutive errors model, is proposed. The optimal prediction in the least squares sense is derived and its properties are studied. Because of its complexity, this estimator needs in practice to be numerically approximated. An automated algorithm is proposed in this thesis. This method is based on an ordering of the observations with respect to the specified spatial process of interest, with the values least in agreement being included towards the end of the ordering. It proves to be useful in case of masked multiple outliers or nonstationary clusters. Simulations are carried out to illustrate its performances and to compare it to other forecasts, robust or not. An application to real data is provided as an illustration of its practical usefulness. The second part of this work also deals with the presence of spatial heterogeneity. One could say that the proposed model offers a characterization of this heterogeneity rather than estimating the locations and sizes of outliers. It is based on the theory of bidimensional α-stable motion. This represents a generalization of the unidimensional Brownian motion. In particular, the stability parameter α can be seen as a measure of the distance between the observations and the hypothesis of a Gaussian distribution. A method of estimation for the parameters of such a process is presented, based on a numerical constrained optimization of the likelihood. Its performances are illustrated by means of simulations. An application ends this second part.