International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Vol XXXV, Part B3. Istanbul 2004 
Their nature is rather different from models usually 
implemented in time-series analysis. This is mainly due to the 
fact that, while the natural flow of time from past to present to 
future imposes a natural ordering or direction on patterns of 
interaction, a two-dimensional model generally possesses no 
such equivalent ordering. 
Hence, our filtering algorithm works under the hypothesis that 
LIDAR measures of terrain/objects height can be rightfully 
represented by means of SAR models. At this step of research, 
first order isotropic SAR models have been employed. 
To introduce SAR models, one can start from the very simple 
expression of a n-dimensional measurement (observation): 
z=p + 8 (1) 
where (specialising it for LIDAR data): 
e zis the [n x 1] vector of LIDAR height values (being n the 
total number of points to be filtered), 
e pn isthe [nx 1] vector of “true” terrain height values, 
e ¢ is the [n x 1] vector of independent and normally 
distributed errors (noise) with mean 0 and variance g2 : 
Considering now for height errors, the effect of a global 
interaction and the spatial dependence among points, height 
values (1) can be rewritten as: 
z=n+(1-pW)'e (2) 
where: 
e Listhe [nx n] identity matrix, 
e pisa value (constant for the whole dataset) that measures 
the mean spatial interaction between neighbouring points, 
e Wisa[nxn] spatial weight (binary) matrix defined as: 
where wj —1 ifthe points are neighbours, wy; =0 otherwise. 
In a regular grid scheme (lattice), a common definition of a 
neighbourhood structure is that wy =1 if the j-th point is 
immediately North, South, East or West of the i-/ point. 
But, since to grid LIDAR data for a lattice scheme leads to a 
loss of information, we preferred operate with raw data, that is 
non-lattice displacement points. In such a non-regular case, W 
can be computed by means of one of the following methods: 
e Distance functions: each LIDAR height measure is linked to 
all others within a specified distance; 
e Nearest neighbor approach: each LIDAR measure is linked 
to its k (k= 1, 2, 3, ...n) nearest neighborhoods; 
e Gabriel’s graph: two generic LIDAR measures i and j are 
said to be contiguous if and only if all other measures are 
outside the i-j circle, that is, the circle whose circumference 
i and j are at opposite points; 
e Delaunay triangulation: all LIDAR measures which share a 
common border in a Dirichlet partitioning of the area are 
joined by an edge. 
This last option was chosen, so obtaining a matrix W with no 
more binary digits but real numbers (Pace, Barry and Sirmans, 
1988); furthermore, a row standardization of W is applied by: 
  
   
   
  
   
   
   
   
  
   
  
  
    
  
  
   
     
   
   
   
    
   
  
  
  
  
  
  
    
    
   
  
  
  
  
  
   
  
  
   
   
  
  
   
   
   
   
  
   
   
   
    
   
  
   
  
   
  
  
  
   
   
n 
W= wi) Zw 
Fl 
1xj 
For (2) to be meaningful, it is assumed that (I— pw)! exists; 
this condition leads to restrain the range of values of p in the 
interval 0 <p < 1. 
An important task is given from the modelling of p, containing 
the “true” height terrain values of (1) by means of some 
polynomial function of East-North coordinates, in such a way to 
analytically define the so-called “trend surface”: 
p=A8 ©) 
where: 
e A isa[n xr] matrix, with A; =f EN 0E) Ni] 
as rows and where s = (r-1)/2, 
e 0= [60 0, 0-1 Tr is a [r x 1] vector of coefficients. 
Equation (3) represents the classical autoregressive problem: 
the estimation of trend coefficients 0 by the measured points. 
Substituting (3) in (2), the general SAR model finally arises: 
z—A0 - (I-oW)'!e (4) 
This form shows the main characteristic of SAR models: they 
require/permit the simultaneous autoregressive estimation of 
both trend 0 and interaction p process parameters. 
Moreover, writing (4) explicitly, we obtain: 
Nj 
Z; =}; +Ej +PL W; (2; -H;) (5) 
j 
with N; number of sites joined to i by an edge (its neighbours). 
Following equation (5), the i-th measured height value can be 
understood as the algebraic sum of two terms: the stochastic 
surface modelling on i-th point (p; 4 £; ) and the global effect 
of errors of such a stochastic modelling on its surrounding 
points (z; —p ; ) via the spatial interaction p. 
3. ESTIMATION OF SAR UNKNOWNS AND OUTLIERS 
Estimation of a spatial autoregressive process with dependent 
variables can be taken over through different approaches. 
A first problem related to this task is due to the computational 
dimension: dealing with millions of LIDAR measures, requires 
great amount of memory storage. For our method, 
computational counts for operation, such as determinants and 
inverses seen in (4), grow with the cube of n. For this 
computational problem, the acquired strip laser data has to be 
suitably shared in sub-zones of about 15.000 points each. 
Moreover, as an analytical problem, traditional maximum 
likelihood techniques require non-linear optimization processes 
using either analytic derivatives or finite difference 
approximations. Unfortunately these can fail to find the global 
optimum and do so without informing the user of their failure. 
Hence summarising, an ideal spatial estimator would: 
e handle large datasets, 
e handle point estimation and inference quickly, 
e not rely on local non-linear optimization algorithms. 
International A 
In order to d 
section, the sc 
estimate unkno 
3.1 Maximum 
For our purpx 
estimate of un 
start from the g 
1(8,p,0%) = (2x 
where: 
is the weigl 
unfortunately 
estimations, he 
as p. It is then 
to 0 and o? ;1 
It can be perf 
by selecting a ' 
and considerin 
L(9, p, 0”) = 
where: 
e eg are the 
regression 
e ey are the 
Thus, to maxir 
L(&.p,o^) 
L(8,p,,0”) | 
L(8,pr, 0”) 
and the value 
among those ir 
The use of a fi 
the chosen va 
make the gram 
estimated py 
precision relat 
log-likelihood 
main property 
Once so obtz 
unknowns and