A MATLAB GEODETIC SOFTWARE FOR PROCESSING AIRBORNE LIDAR BATHYMETRY DATA

The ability to build three-dimensional models through technologies based on satellite navigation systems GNSS and the continuous development of new sensors, as Airborne Laser Scanning Hydrography (ALH), data acquisition methods and 3D multi-resolution representations, have contributed significantly to the digital 3D documentation, mapping, preservation and representation of landscapes and heritage as well as to the growth of research in this fields. However, GNSS systems led to the use of the ellipsoidal height; to transform this height in orthometric is necessary to know a geoid undulation model. The latest and most accurate global geoid undulation model, available worldwide, is EGM2008 which has been publicly released by the U.S. National Geospatial-Intelligence Agency (NGA) EGM Development Team. Therefore, given the availability and accuracy of this geoid model, we can use it in geomatics applications that require the conversion of heights. Using this model, to correct the elevation of a point does not coincide with any node must interpolate elevation information of adjacent nodes. The purpose of this paper is produce a Matlab® geodetic software for processing airborne LIDAR bathymetry data. In particular we want to focus on the point clouds in ASPRS LAS format and convert the ellipsoidal height in orthometric. The algorithm, valid on the whole globe and operative for all UTM zones, allows the conversion of ellipsoidal heights using the EGM2008 model. Of this model we analyse the slopes which occur, in some critical areas, between the nodes of the undulations grid; we will focus our attention on the marine areas verifying the impact that the slopes have in the calculation of the orthometric height and, consequently, in the accuracy of the in the 3-D point clouds. This experiment will be carried out by analysing a LAS APRS file containing topographic and bathymetric data collected with LIDAR systems along the coasts of Oregon and Washington (USA).


INTRODUCTION
A direct georeferencing (DG) system provides the ability to directly relate the data collected by a remote sensing system to the Earth, by accurately measuring the geographic position and orientation of the sensor without the use of traditional groundbased measurements (Mostafa et al., 2001;Pepe et al., 2015).Examples of where DG systems are used in the airborne mapping industry include: Airborne Laser Scanning (ALS) (Mallet et al., 2009) and Airborne Laser Scanning Hydrography or LIDAR bathymetry (ALB) (Mandlburger et al., 2013).Because the technique is based on DG ellipsoidal heights, the coordinates must be converted to orthometric system through a suitable geoid undulation model.The grid structure is necessary because the geoid undulation N (, λ), being it a function of the position, cannot be expressed in analytical form closed but only by numerical approximations.The equation relates the orthometric height H (above the geoid), the ellipsoidal height h and the geoidal undulation N, is (Hofmann et al., 2006): The undulation value is computed by interpolating over the corresponding geoid model grid.Studies conducted on Alpine region on the influence of interpolation techniques in correcting ellipsoidal in orthometric height have showed that the "Triangulation with linear interpolation" (Pepe, 2014) and "bilinear" (Pérez et al., 2012), allow you to process millions of points in a few seconds and without changing the value of the interpolated point (Pepe, 2014).With recent advances in LIDAR (Light Detection And Ranging) technology, coastal managers now have the capacity to acquire high resolution digital elevation data covering the littoral or inter-tidal zone.This development of hardware has been followed by a development of a new data file format standard known as the American society for Photogrammetry and Remote Sensing (ASPRS) Lidar Exchange Format (LAS) (Samberg, 2007).The ASPRS LAS is a binary file format public, for the exchange of three-dimensional point clouds, which stores information about the types of data; it, despite it was developed primarily for the exchange of data , point clouds LIDAR, has also become a standard for the exchange of any file in the format x, y, z.This format, created as an alternative to proprietary files produced by the manufacturing companies of the LIDAR and to those who use file-type ASCII, solves two great problems.The first is related to performance, because the reading and interpretation of the elevation data ASCII can be very slow and the resulting file size can be extremely large, even for small amounts of data.Empirical tests, carried out on ASPRS LAS files (LAS 1.1) that contain information such as coordinates, intensity, return number of returns (X,Y,Z,I,R,N), have shown that the size of these files are reduced by about 55% compared to a simple ASCII file (Pepe, 2014) (Quintero, 2013).The second problem is to maintain all the specific information of the LIDAR data; this is particularly useful in applications of airborne laser scanning as in them, working on large areas, the geoid undulation can take considerably variable values in the survey area; transform the point clouds with a height orthometric preserves information about the nature and quality of the point.However, in the literature there are software solutions (Hug et al., 2004) but who do not are not in the public domain and not allow you to check for errors or implementations as may occur with this algorithm in Matlab®.

EGM2008 EARTH GEOPOTENTIAL MODEL
The latest and most accurate global geoid undulation model, available worldwide, is EGM2008 (Pavlis et al., 2008).From NGA website it is possible to download the global grid of the geoid undulation values, in geographic coordinates, with a resolution of 2.5'x2.5'.This model is available online at the website of the NGA (Figure 1): http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/ Figure 1.Global 2.5' x 2.5' geoid undulations This grid was formed by merging terrestrial, altimetry-derived, and airborne gravity data.Over areas where only lower resolution gravity data were available, their spectral content was supplemented with gravitational information implied by the topography.EGM2008 is complete to degree and order 2159, and contains additional coefficients up to degree 2190 and order 2159.Over areas covered with high quality gravity data, the discrepancies between EGM2008 geoid undulations and independent GPS/Levelling values are on the order of ±5 to ±10 cm (Pavlis et al., 2012).In the Newton's Bulletin ( 2009) several studies on the accuracy and quality of the EGM2008 model in different parts of the world are described.In general, the different models of geoid undulation are provided in the form of grids, ie files containing a table of n values calculated on the nodes of a predetermined grid.Many of the studies in the literature about the steep slope of the undulations of the geoid which occur on mountain ranges, but strong values of slopes are realized even outside of the land emerged; just consider that along the coasts of Indian Ocean there are slope values of about 0.15‰.This slope value, for example, for a length 30 km will cause a variation of the geoid undulation of 4.5m.
A profile of the geoid undulation close to the Sri Lanka coast is shown in Figure 2.
Figure 2. Profile of the geoid undulation close to the Sri Lanka coast

STRUCTURE OF ALGORITHM
The specific software, for processing millions of points, divide the point clouds obtained by aerial survey in n geographic blocks; subsequently, with special routines, they classify the ground from the other elements (buildings, trees, etc.).From this procedure, you can get different cartographic products, such as DTM, DSM or DEM, in which reference is made to ellipsoidal heights; is, therefore, necessary to convert ellipsoidal heights to orthometric.Usually the orthometric height is determined from the n blocks by performing n conversions by using other specific software.For example in Italy to carry out the transformation of heights it is used Verto software, which allows you to process files whose structure is necessarily of the type: north, east, ellipsoidal height.Instead, if we convert the heights ASPRS LAS file directly, we will use a single height conversion operation.From here the idea of implementing, in Matlab® environment, a particular algorithm that allows, according to the diagram shown in Figure 3, the conversion of heights.

Figure 3. Flow chart of algorithm
The basic structure for reading files ASPRS is due to Alexander (2008).The input file is modified by the function "change z" which changes the value of the height and because this function change the basic structure of the ASPRS LAS file it is necessary to make a copy of the original file before you call the function.
For the transformation of heights, the bilinear interpolation is used; this interpolation considers the four grid values surrounding the point to be interpolated.The interpolated final point is computed with two successive linear interpolations, for line and column, using the undulation values of the two surrounding points the point interpolated in each direction.This function is implemented in MATLAB by means the module "interp2".
The algorithm, valid on the whole globe and operative for all UTM zones, allows the conversion of ellipsoidal heights using the EGM2008 model.If the coordinate of point clouds are in geographic coordinates, in the software there is an algorithm that allows to transform the by geographic coordinates to plane coordinates, through the well-transform known in the literature Having available a more accurate geoid model, it can be used provided it is in the format East, North, geoid undulation.

CASE STUDY
This experiment will be carried out by analysing a LAS APRS file containing topographic and bathymetric data collected with LIDAR systems along the coasts of Oregon and Washington (USA) from the NOAA Coastal Services Center (http://www.csc.noaa.goc/crs/tcm/index.htm)for our study.This data set published by the NOAA's Ocean Service, Office for Coastal Management (OCM), is an LAZ (compressed LAS) format file containing LIDAR point clouds data.These files contain topographic and bathymetric LIDAR data collected with Leica ALS60 (topography) and SHOALS-1000T (bathymetry) systems.
To transform the ellipsoidal to orthometric heights, with reference to the North American Vertical Datum of 1988 (NAVD88), the GEOID09 model was used.The position data were obtained using post processing KGPS methods.Because this data set includes several Bathymetry and Topography, the vertical accuracy of the data is better than ±15 cm RMSE.The vertical accuracy of topographic LIDAR is dependent on several external factors (Turton, 2006); to account fully for the influence of all potential error sources, the 1 sigma vertical accuracy of topographic LIDAR is generally quoted to be ±0.15m(Quadros et al., 2008).
In the Figure 4 is shown an RGB image of the study area which covers part of the Oregon coasts and more precisely: 44°15'54" ≤ φ ≤ 44°18'54" 124°06'00" ≤ λ ≤ 124°07'00" An example of these bathymetric LIDAR data encoded with LAS 1.2., showed using the software Global Mapper v.15.0 (Figure 5) from which it can be noted that the elevations of about 276.000 points ranging from -20m to 35m.The structure of this ASPRS LAS file is reported in table 1.

FIELD VALUE
x  The transformation of the height (elevation) of this point clouds was performed in about 1 second.

Figure 5. Example of bathymetric LIDAR data
The geoidal undulation value varies, along the strip, from -23.69m and -23.01m, with a difference of 68 cm.Consequently, using a constant value to correct the strip height you make an error higher than the accuracy obtained with an airborne LIDAR.
Of course the LAS file, from the NOAA Coastal Services Center, already contains the orthometric height.This operation was conducted in order to verify the effective functioning of the algorithm on a sample of bathymetric LIDAR.

CONCLUSIONS
The software developed in Matlab® allowed to transform the point clouds from ellipsoidal to orthometric height.The value of the geoid undulation was taken from the geoid model EGM2008, thanks to its high performance in terms of accuracy.In relation to the precision obtainable with LIDAR systems, theoretically, when the gradient of the undulations does not exceed 15 cm you may use a constant value of the undulations; this means that in areas with steep slopes of geoid undulation, it is suggested to shorten the length of the strip, of the survey with airborne LIDAR, in order not to exceed the threshold of accuracy required.However, given the speed of processing we always suggest the use of interpolation for the correction of the height.Moreover we hope a greater use of the LAS format.1.4 because still little used.In fact, the most recent US Army Corps of Engineer data (available online) was encoded with LAS 1.2.Lastly, the software developed is useful, not only in bathymetric applications, but in all geomatic applications that produce point clouds which require the corrections of the ellipsoidal height in orthometric.

Figure
Figure 4. Study area [276542x1 double] y [276542x1 double] z [276542x1 double] Intensity [276542x1 double] ReturnNumber [276542x1 double] NumberOfReturns [276542x1 double] Classification [276542x1 double] ScanAngleRank [276542x1 double] ScanAngleRankInfo '-90 to 90 left side' UserID [276542x1 double] PointSourceID [276542x1 double] TimeInfo 'Time is "standard GPS time minus 1e9"' Time [276542x1 double] . This is particularly useful especially in file structures containing multiple fields such as the one in version LAS 1.4.In early August, 2013, the LIDAR Division introduced the concept of a LAS Domain Profile, enabling expansion of the base LAS 1.4, such that specialized communities can extend the format to their needs.The first approved and published LAS Domain Profile was designed to support the coastal LIDAR community.The Topo-Bathy LIDAR Domain Profile adds point classification values for bathymetric point, water surface, derived water surface, submerged object, IHO S-57 object, and bottom not found depth.Extra Byte Variable Length Records (EXTRA_BYTES or Extra Byte VLRs) are added for pseudoreflectance, uncertainty, water column depth, Figure of merit, and processing specific flags

Table 1 .
Structure of this ASPRS LAS file