variable at unsampled sites within the area covered by sample 
locations. Normally, three kinds of approximation are often 
used. (1) If there is only one point in one grid, then the point is 
used to represent the value of the grid. (2) If more than one 
point in one grid, the average value of the points is used, or the 
median, maximum, or minimum one is chosen. (3) If no data 
point is in the grid, the nearest point to the grid centre is used. 
In conclusion, the above methods all use the value of a certain 
point to represent the grid of interest. This scheme usually leads 
to part of the natural dispersion not reflected by the data, and 
causes information loss. In order to better represent the data, 
geostatistic approach for spatial interpolation is considered by 
modelling the spatial correlation among the data points in a 
certain neighbourhood. 
A fundamental assumption for geostatistical methods is that any 
two locations that are a similar distance and direction from each 
other should have a similar difference squared. This relationship 
is called stationary. If the spatial process is isotropy, spatial 
autocorrelation may depend only on the distance between two 
locations. The rate at which the correlation decays can be 
expressed as a function of distance. If the process is second 
order stationary, the covariance between any two random errors 
depends only on the distance and direction that separates them, 
not their exact locations. Semivariogram is a common tool to 
capture the second-moment structure of spatial data. It describes 
the variability of two locations of the data separated by a 
distance h (Oliver et al., 2005): 
y(s,h) = -Var[Z(s)-Z(s + hj\ (1) 
Where y(s ,h) is the semivariogram, and Z(s) is the data value at 
location s. 
Wallace and Marsh (2005) used geostatistics to extract 
measures that characterize the spatial structure of vegetated 
landscapes from satellite image for mapping endangered 
Sonoran pronghorn habitat. They use variogram parameters to 
discriminate between different species-specific vegetation 
associations. Woolard and Colby (2002) used DEM generated 
from ALS data and spatial statistics to better understand dune 
characterization at a series of spatial resolutions. 
Digital images are rich in data, but in many instances they are 
so complex as to require spatial filtering to distinguish the 
structures in them and facilitate interpretation. The filtering can 
be done geostatistically by Kriging analysis. It proceeds in two 
stages. The first involves modeling the correlation structure in 
an image by decomposing the variogram into independent 
spatial components. The second takes each component in turn 
and kriges it, thereby filtering it from the others. 
In Lloyd and Atkinson (2002), inverse distance weighting, 
ordinary Kriging and Kriging with a trend model are assessed 
for the construction of DSMs from ALS data. Factorial Kriging 
is a geostatistical technique that allows the filtering of spatial 
components identified from nested variograms. 
random field, and denotes by the data vector, Z(s) = {Z(si), ..., 
Raster representation of the random field means to lattice the 
continuous domain B and then calculate a typical elevation for 
each grid (pixel). Let the region of a given pixel be B and the 
corresponding area be |5|. The elevation of the pixel can be 
predicted by calculating the average value of the random field 
in B: 
Z{B) = ~Y^Z{s)ds 
(2) 
Since the exact value of Z(B) cannot be calculate directly, we 
can only predict Z(B) using the observed data. Therefore, a 
window is set up around the grid B, raw data points 
s* = K e s;s k e B ) 
(3) 
Then, the window B are used to predict Z(B). With block 
Kriging predictor: 
/>(Z,Z(5)) = £V^) 
(4) 
where A, k are chosen to minimize the mean-squared prediction 
error between p(Z,Z(B)) and real, unknown elevation of B: 
e[(e[z(B)\Z(s)]-Z(B)J\ 
(5) 
Under this circumstance, p(Z,Z(B)) is an unbiased prediction of 
Z(B), the optimal weights {X k } can be obtained by: 
(<7(1,■»))+!— 
(6) 
where the elements of the vector c(5, s) are Cov[Z(5), Z(s*)]. 
By discrentizing B into points, {pj}, the point to block 
covariance can be approximated using 
Cov[Z(B),Z(s)\ * 1/N^du'j,s) 
Where £ is the matrix composed by the covariance of every 
two observed points. Due to the autocorrelation of spatial 
process, the elevation values in a small region always have a 
constant mean u(s), the covariance between observations Z(s,), 
Z(sj) is: 
Cov(Z(s i ),Z(s J )) 
= E[{Z(s,)-u(s)}{Z(s J )-u(s)}] 
As shown in figure 2, block Kriging is a method which uses the 
value of a block to represent the value of the grid.