**Nishank Agrawal, Nitish Kumar Agrawal, Bharat Lohani**

Department of Civil Engineering

Indian Institute of Technology Kanpur, Kanpur 208016 (India)

Tel: 259 7623, Fax: 259 7395

Email: [email protected]

**1. Introduction:**

Laser altimetry is a direct method to efficiently map the earth's surface in three spatial dimensions. Topographic mapping with LiDAR seems on the verge of supplanting traditional photogrammetry in the near future for detailed terrain mapping. This technique has a large and growing number of actual and anticipated applications in disciplines ranging from urban analysis to natural resources. Initially the focus of research in this field was on how to obtain terrain geography and objects such as buildings, trees etc. Furthermore, several other researchers were interested in identifying the error involved in LiDAR data. In all these studies raw LiDAR data were employed and attempts were made to compare results with ground truth. For a detailed review of these efforts please refer Lohani (1999) and special LiDAR issue of ISPRS Journal (ISPRS, 1999). However, the literature clearly indicates that a limited success has been achieved in these efforts and a lot needs to be done to make full utilization of this technology. The basic drawback of methods relying only on field data is that, it is not always possible to collect true ground truth for validation of these approaches. In view of the above it is of utmost importance that research be carried out to develop theoretical equivalents of LiDAR data collection process. Realizing this, recently the research related to simulation of LiDAR data has started taking place. The purpose of these researches should be to study the effect of various flight parameters (flying height, direction of flight etc.), LiDAR instrument parameters (firing frequency, scanning angle, number of scan lines), wrong instrument calibration and environment on LiDAR data generated for different types of objects and terrains being scanned. Beinat and Crosilla (2002) created a virtual scenery with four distinct buildings using MatlabTM and then produced some LiDAR images by randomly sampling the surface points of the model with an average resolution of 0.5 points/m2 and a random error of ±10 cm. In this way the effect of different weighting models was studied. Holmgren et al. (2003) described the effect of LiDAR scanning angle for estimation of mean tree height and canopy closure. For this a computer modeled forest was simulated for different scanning angles. Most of these researchers have studied the effect of one parameter on one object only, and it was assumed that either the effect of other parameters is negligible or independent of the effect of parameters being studied. However, in real world all parameters affect together and have varying degree of effect on different terrain conditions. Furthermore, as in case of Beinat and Crosilla (2002) the data points were assumed to be randomly distributed which does not reflect the way LiDAR collects data.

**2. Outline of present work:**

The aim of present work is to develop a simulator which will be helpful in simulation of more than one parameter on surfaces which will resemble terrain in real. A simulator is a computer programme that creates the model of existing external conditions and generates output. In the present case external conditions are LiDAR data collection system, and output generation is data points. The simulation is realized using a set of mathematical equations and coding them in C programming language. With the help of simulator at any instant the location of aircraft and laser vector angles with three axes are obtained, thus equation of laser vector is known. Surfaces with known mathematical equations are chosen to represent ground. More than one mathematical surfaces are stitched together to resemble complex terrain. Data points are obtained by the intersection of laser vector and the surface (Figure 1). Nature of surfaces is chosen with the aim to model the possible ground types and the objects such as buildings, trees etc. The data points thus obtained can be plotted in suitable software (Matlab here) to visualize the output of simulator. Various such displays can be obtained to study the effect of various parameters on different terrain surfaces.

**Figure 1: Working of LIDAR simulator **

The specifications of ALTM 3025 by Optech Inc. have been used in the present case, which gathers up to 25000 points per second, with maximum flying altitude of 3000 m., swath angle varies from 0° to ±20°, and the number of scans from 0 to 70 in one second (ALTM 3025 specifications, . However, with a little change other sensors and specifications can also be simulated.

**3. Methodology for simulation**

The simulation in case of LiDAR can be categorized in two distinct parts 1) sensor performance and environment 2) Nature of terrain. The former in the present case is implemented by considering a straight line path of aircraft between two points, input by user, for starting and finishing position of aircraft which in real are given by a kinematic differential GPS, mounted on top of aircraft. User also defines the swath angle and number of scans desired to carry out simulation. In accordance with the above defined parameters, the instantaneous coordinates of aircraft and direction of laser beam are calculated. Equation of laser vector thus is given by:

Where, and ω, γ and κ are absolute angles with positive x, y and z axes respectively, and X0, Y0, Z0 are instantaneous coordinates of aircraft at the time of firing the pulse (Figure 1). In the present case it is assumed that there are no roll, pitch, and heading movement of aircraft, though the same will be accounted for in future.

The second part of simulation is the nature of terrain, for which LiDAR data is being gathered. This terrain can be represented by mathematical surfaces. The following paragraphs describe the procedure for generating LiDAR data on these surfaces. Surfaces ranging from simple planar surfaces to complex i.e. combination of multiple surfaces have been discussed. The simple surfaces have equations of the form of having no border limitations. X and Y from equation (2.1) are put into the equation of surface and a simple equation of the form of is obtained which is solved for Z. Now from equation (2.1) again the value of X and Y corresponding to Z obtained above is calculated i. e.

The surfaces of the above form are:

1. Plane surface (Ax +By +Cz +D =0) (Figure 2.1),

2. Cone (z2 =x2+y2) (Figure 2.2),

3. z =x2+y2 (Figure 2.3).

After this, to involve a little more complexity, two or more surfaces are stitched together:

- 1. Stitching of two planes (Figure 2.4).

The intersection of laser vector with both planes is taken and the higher value of Z is chosen and data points from corresponding plane are obtained.2. Inverted cone on horizontal surface (Figure 2.5).

The intersection with cone is obtained similar to previously discussed cone except that, here altitude of the cone top is mentioned by user along with the reduced level of surrounding ground. If Z value obtained with the intersection of cone is less than that of ground reduced level, then intersection with the ground is sought.3. Hemisphere ( ) on horizontal ground (Figure 2.6).

If D<0; it indicates no intersection with sphere and thus the intersection with

ground surface is sought.

else, , where R is the radius of hemisphere.4. Right Circular Cylinder on horizontal surface:

If D<0: no intersection with cylinder, therefore Z is obtained by intersection with ground.

else

- If Z&glt; Height of the cylinder, then actual point is either on top of cylinder or on ground, put Z= Height of the cylinder, obtain X and Y.
If then point is on ground (obtain Z by intersection with ground) else Z= Height of the cylinder (On top of cylinder).

If Z>ground RL and < Height of the cylinder, then point is either correct or at top of cylinder, put Z= Height of the cylinder. Obtain X and Y.

- If then Z obtained before else Z= Height of the cylinder.
If Z

- If then point is on ground (obtain Z by intersection with ground) else Z= Height of the cylinder (On top of cylinder).
5. Inverted cone/ hemisphere on top of cylinder placed on horizontal surface

Intersection with cone/ hemisphere (Figure 2.7 and 2.8) is obtained similar to as explained before, but when Z obtained is less than or equal to cylinder height, the laser vector is solved with cylinder as discussed above.

To represent variability of natural terrain more complex surfaces are chosen. These surfaces have an equation of the form of . To solve the laser vector and surface, numerical methods are employed since no exact mathematical expression can be obtained in this case for solving X, Y or Z.

Here Z and Y are replaced in terms of X from equation (2.1) in surface equation and an equation of the form of is obtained. Fixed point iteration method is used here. Simulator chooses starting point of iteration such that , which indicates that the iteration will converge (Chapra and Canale, 2000), an accuracy of 0.1 μm is obtained before terminating iteration.

Surfaces chosen of this form are:

**Figure 2: Surfaces generated by simulator**

**4. Error Modeling:**

After modeling for above ground and instrument properties, random errors in vertical as well as horizontal directions are included in simulated data. It has been observed that these random errors follow a normal distribution with a standard deviation of around 15 cm and 50 cm in vertical and horizontal directions respectively (Lohani, 1999).To realize this, the rand() function of C programming language is used (Gottfried, 1988). Before calling this function another function srand(time(NULL)) is used, which sets its argument as the seed for a new sequence of random numbers to be returned by rand(). These functions are included in the library stdlib.h. Random numbers following normal distribution are generated by the procedure as described above with appropriate standard deviation (s) and stored in an array of length 25000 which is same as the number of points being scanned per second. Then for each point generated by the intersection of laser beam and surface these values are added, so that every point has unknown amount of error, independent of each other. This kind of error can not be removed at this point as source and amount of these are unknown for every point, necessitating special measures for their minimization.

**5. Results and discussion:**

The results in following paragraphs have been presented for only two kinds of surfaces in view of the space limitation. The points obtained as a result of the intersection of the laser vector and the surfaces are plotted using MatlabTM. This surface is then visually compared with the original surface. The comparison at present aims at ensuring accuracy of simulator as data points simulated will form the same geometry as of original surface. Furthermore, graphical results are generated to show the effect of relative positioning of aircraft flight line and ground on the shadow.*5.1 Example 1.*

An attempt has been made to represent a city model by placing cylinders of different radii and heights at arbitrary locations on a plane surface, the number, radii and heights of the cylinders being user defined. The aircraft is then flown as per user's choice of altitude, flying direction, scanning frequency, firing frequency.For following user defined parameters (Table 1), the display was obtained:

Flight in one second from (-200, 0, 1500) to (200, 0, 1500),

swath angle 20° from mean,

number of scans = 50,

neighboring ground with reduced level = 0 m.

**Table 1: Cylinder data for model 1**No. Center Radius Height 1 -60, -10 20 300 2 -30, 70 20 350 3 -40,-200 30 250 4 -5, 400 15 185 5 20, -150 20 225 6 55, 300 20 325 7 60, -400 15 300 Following display was obtained (Figure 3):

**Figure 3: City model with different cylinders**

As can be seen here, there are different densities of data points on curved surfaces of cylinders. There is more data point density on the curved surfaces of the cylinders nearer to the flight path. It can also be observed that there are no points on the rear side of the cylinders as the laser beam is unable to reach there. There have been no data points on the curved surface of cylinder '1' as the aircraft is flown almost exactly above this cylinder as per the given flight parameters in the present case.When a laser beam is returned back from an object on the terrain, it can not reach some portion of terrain adjacent to that object lying on the leeward side. Hence no data points are obtained in that region. This is similar to the shadow of that object, and is called Shadow Zone.

It can be observed that the shadow zone corresponding to the cylinders increases as the height of the cylinders and their distance from the flight line path increases.

*5.2 Example 2.*

For studying the effect of distance of object from flight line path on shadow zone, the following model was generated. Data for surface are given in Table 2.

**Table 2: Cylinder data for model 2**Number Center Radius Height 1 -50, -50 30 30 2 -50, 50 30 30 3 50, -50 30 30 4 50, 50 30 30 Flight line is varied in three cases, as described in figure 4.

**Figure 4: Variation in shadow zone with variation in flight lines.**

Case 1: Flight direction along X axis, Y=0 Case2: Flight direction along X axis, Y=25 Case3: Flight direction along X axis, Y=150Figure 4: Variation in shadow zone with variation in flight lines.When aircraft is flying above X axis (case1), it is symmetrical to four cylinder model and shadow zone of each of four cylinder is observed to be identical. When flight line is moved to a different location parallel to X axis (case 2), the shadow zone of those cylinders which are at more distance is enlarged as compared to those which are nearer. In case 3, when the flight line is moved further away from X axis, the shadow zone of farther cylinders has enlarged to a more extent, and the shadow zone of nearer cylinders has shifted its direction as laser beams are no more obstructed by object to reach previous shadow zone.

**6. Future course of work:**

This paper describes the progress achieved so far in the design of simulator and accordingly a selected few results are presented. It is planned to carry out the following in rest of the project:- Inclusion of more sources of error, systematic and random, resulting due to sensor functioning and malfunctioning.
- Inclusion of error due to aircraft trajectory perturbations, measured by Inertial Navigation System (INS) (Shukla, 2003).
- The results obtained by simulator will be evaluated quantitatively.

Furthermore, it is required to study the data simulated for understanding the effect of change of various parameters.

**7. Conclusion:**

The simulator described in this paper, is a user friendly programme, where user can define various flight parameters and terrain conditions as per his scope of survey. The data points obtained from the simulator can be drawn in suitable 3D plotting software for varying flight/ instrument parameters (more than one at a time) with different terrain conditions, and thus their effect can be studied.The main application of this simulator would be in economizing the LiDAR data collection process. Before actually flying the aircraft for data collection, the terrain condition resembling to actual ones, can be modeled in laboratory, with the help of this simulator for different flight parameters. While collecting the data for a dense forest less density of data points is sufficient, but when data for city models is collected high data point density is required to efficiently map the features. A particular object and its neighboring area may be more important for the purpose of data collection, thus might require lesser shadow zone. Therefore, several outputs of simulator with different parameters can be compared to obtain a particular set of parameters for which required information can be obtained in the most economical way. Hence the decisions regarding the flight parameters would no more be taken on the basis of thumb rule or past experience of the surveyor, but on facts based on solid scientific thinking.

This simulator can also be used to assess the accuracy of existing information extraction algorithms. For this, the existing algorithms can be operated on simulated data and information extracted can be compared with the true mathematical surfaces forming the terrain.Furthermore, the simulator is to accommodate different kinds of systematic and random errors, thus a better understanding of these errors with different instrument and environmental conditions can be achieved and this may help in further development of error removal techniques.

Lastly, the simulator can also play a significant role in education and training on LiDAR technology.

**References:**- ALTM 3025 specifications, URL:
- Axelsson, P., 1999. Processing of laser scanner data- algorithms and applications, ISPRS Journal of Photogrammetry & Remote Sensing 54 (1999) 138-147
- Beinat, A., Crosilla, F., 2002. A Generalized factored stochastic model for the optimal global registration of LiDAR range images,
- Chapra, S. C., Canale, R. P., 2000. Numerical Methods for Engineers With Programming and Software Applications, Tata McGraw-Hill Publishing Company Limited, New Delhi, India, 145 p.
- Gottfried, B. S., 1988. Programming with C, Tata McGraw-Hill Publishing Company Limited, New Delhi, India, 183 p.
- Holmgren, J., Nilsson, M., Olsson, H., 2003. Simulating the effect of lidar scanning angle for estimation of mean tree height and canopy closure, Canadian Journal of Remote Sensing, Vol 29, No. 5, pp. 623-632, 2003.
- ISPRS, 1999. ISPRS Journal of Photogrammetry and Remote Sensing volume 54, No 2-3
- Lohani, B., 1999. Application of Airborne Remote Sensing to the study of Intertidal zone, Ph. D. thesis, The university of Reading, Reading, UK
- Maas, H. G., Vosselman, G., 1999. Two algorithms for extracting building models from raw laser altimetry data, ISPRS Journal of Photogrammetry & Remote Sensing 54 (1999) 153-167
- Shukla, R. K., 2003. Inertial Navigation Systems, Short term course on Airborne Altimetric LiDAR, IIT Kanpur, 10-14 March 2003.

- If Z&glt; Height of the cylinder, then actual point is either on top of cylinder or on ground, put Z= Height of the cylinder, obtain X and Y.