High spatial resolution digital elevation model ( DEM ) production using different interpolations techniques

DEMs, thus, simply regular grids of elevation measurements over the land surface.The aim of the present work is to produce high resolution DEM for certain investigated region (i.e. Baghdad University Campus\ college of science). The easting and northing of 90 locations, including the ground-base and buildings of the studied area, have been obtained by field survey using global positioning system (GPS). The image of the investigated area has been extracted from Quick-Bird satellite sensor (with spatial resolution of 0.6 m). It has been geo-referenced and rectified  using 1st order polynomial transformation. many interpolation methods have been used to estimate the elevation such as ordinary Kriging, inverse distance weighted (IDW) and  natural neighbor methods. The mosaic  algorithm has then been applied between the base and building layers of studied area in order to perform the final DEM. The accuracy assessments of the interpolation methods have been calculated using the root-mean-square-error (RMSE) criterion. Finally, the estimated DEMs have been used to constructing 3-D views of the original image.


Introduction
A digital elevation model (DEM) is a digital representation of ground surface topography or terrain.It is also widely known as a digital terrain model (DTM).A DEM can be represented as a raster (a grid of squares) or as a triangular irregular network.DEMs are commonly built using remote sensing techniques, but they may also be built from land surveying.DEMs are used often in geographic information systems, and are the most common basis for digitally-produced relief maps.
A DEM is used as a mean of 3-D terrain modeling; it serves as a basic source of information for deriving geo-spatial uniqueness.Currently, Digital elevation modeling (DEM) is one of the modern methods for representing the topographic surface of the terrain, (i.e., how the elevation of the ground surface is changing with position).Traditionally this has been done by contour lines on topographic maps.DEMs being generated by many methods, such as ground survey, photogrammetry, Light Detection and Ranging (LIDAR), and Global positioning system (GPS) techniques which have been adopted to generate high resolution digital elevation model [1,2].Where GPS, well known as a versatile, global tool for positioning, has also become the primary system for distributing time and frequency.GPS receivers are fixtures in telecommunication networks, and calibration and testing laboratories.They make it possible to synchronize clocks and calibrate and control oscillators in any facility that can place an antenna outdoors for line-of-sight reception of the GPS satellites.The GPS satellites are controlled and operated by the United States Department of Defense (USDOD).The constellation includes at least 24 satellites that orbit the earth at a height of 20,200 km in six fixed planes inclined 55 from the equator.The orbital period is 11 h 58 m, which means that a satellite will orbit the earth twice per day.By processing signals received from the satellites [3,4].The GPS satellites broadcast on two carrier frequencies: L1 at 1575.42 MHz, and L2 at 1227.6 MHz.Each satellite broadcasts a spread spectrum waveform, called a pseudorandom noise (PRN) code on L1 and L2, and each satellite is identified by the PRN code it transmits.There are two types of PRN codes.The first type is a coarse acquisition (C/A) code with a chip rate of 1023 chips per millisecond.The second type is a precision (P) code with a chip rate of 10230 chips per millisecond.The C/A code is broadcast on L1, and the P code is broadcast on both L1 and L2 [5].GPS reception is line-of-sight, which means that the antenna must have a clear view of the sky.If a clear sky view is available, the signals can be received nearly anywhere on earth.the absolute WGS84 coordinates of at least one site have to be known accurately as all measurements in GPS system are depicted in WGS84 coordinate system.And it consists of three segments space segment see Fig. 1, control segment and user segment [6,7].

Studied area
study area located in the middle of Iraq, Baghdad City center, it represent Baghdad university area, which consists of different buildings include the colleges and its department, vegetation and some of trees cover the land of this region, in this research college of science has been used to built it's digital elevation model and the general information of the original image can be showed in Table 1 and Fig. 2 shows the studied area.

Methodology
this research is depended on splitting the region into two levels; i.e. building and ground levels.The coordinates of the two levels were measured by the ordinary GPS, ninety (90) points have been measured on the ground level of the studied region.Many interpolation techniques are implemented on just the ground elevation points.Reclassify the building layer to specify the real height of each class, after that mosaic model has been used to project the building layer on the ground layer in order to get the final shape of the digital elevation model (DEM), root mean square error has been employed to show appropriate interpolation method, the applying interpolation methods can be explained as follow:

Interpolation Methods
Interpolation methods can be classified into deterministic and geostatistical.Deterministic interpolation methods use mathematical functions to calculate the values at unknown locations based either on the degree of similarity (e.g.IDW) in relation with neighboring data points.Geostatistical techniques use both, mathematical and statistical methods to predict values at unknown locations and to and to provide probabilistic estimates of the quality of the interpolation based on the spatial autocorrelation among data points, such as kriging methods [8,9].

Inverse Distance Weighting (IDW)
Inverse distance weighting models work on the premise that observations further away should have their contributions diminished according to how far away they are.The simplest model involves dividing each of the observations by the distance it is from the target point raised to a power α : The value k j in this expression is an adjustment to ensure that the weights add up to 1.If the parameter α=1 we have: Many GIS packages provide this kind of inverse distance model for interpolation, as it is simple to implement and to understand.Often the model is generalized in a number of ways: a faster rate of distance decay may be provided, by including a power function of distance, α>1, rather than simple linear distance.While any α value convenient for a given application may be used, common practice is to use distance (α=1) or distance squared (α=2).
since it is possible that the grid intersection could coincide with a data point (especially likely at region corners or edges if MBRs have been used on the original point set), an explicit check or adjustment to the expression is needed to avoid computational errors (overflow).Typically the adjacent point weight is set to 1 (its value is copied) and the remaining weights are set to 0. if a user-selectable adjustment to the minimum inter-point distance is specified, this can result in smoothed rather than exact interpolation.This may be a simple incremental amount added to the distance, t, or an adjusted distance value such as: as with all methods, additional controls may be applied or available: limiting the number of points included; specifying the search directions and search shape; and limiting computations by excluding pre-defined regions, break lines or faults.Fig. 2 shows how simple IDW with no smoothing and power 2 distance decay results in dips and peaks around the data points but is otherwise relatively smooth in appearance [10].

Ordinary Kriging:
Ordinary kriging is the most general and widely used of the kriging methods, it is divided into two distinct tasks:  Quantifying the spatial structure of the data (known as variography) and producing a prediction i.e., fitting a spatial dependence model to the data. Make a prediction for the unknown value of a specific location.Achieved by using the fitted model from the variography (spatial data configuration) and values of the measured sample points around the prediction location.The equation used in Ordinary Kriging is [11,12]:-Z* (u) is the Ordinary Kriging estimate at spatial location u, n (u) is the number of the data used at the known locations given a neighborhood.Z (u α ) are the n measured data at locations u α located close to u m= mean of distribution λ α (u)= weights for location u α computed from the spatial covariance matrix based on the spatial continuity (variogram) model, which is given by: n is the number of data pairs separated by distance h z(u i ) and z(u i +h) are the data values at locations separated by distance h.

Geometric Correction (Rectification):
Different polynomials may be used to convert the source image's coordinates to the rectified map's coordinates.Number of the utilized GCPs used depends on the distortion in the imagery.Since the covered area in this paper (i.e.Baghdad University Compass\ college of science) is not too large (can be covered by small size satellite image) for more details see Table 1 , thus the 1 st order polynomial (affine transform) was found to be too efficient for the rectification process.To ensure the coincident of the rectified and resampled image points with the GCPs (as measured by the GPS), we have distributed them on the image plane, illustrated in Fig. 3 and 4.
General, remotely sensed images are gathered by a satellite or aircraft represent the irregular surface of the Earth.Even images of seemingly flat area are distorted by both the curvature of the Earth surface and the sensor being used.In what follows, the mathematical operations concerning the geometrical correction, the goal of geometric correction can be abbreviated by: holding the mismatch to a displacement of no more than one pixel.Under ideal conditions, and with as many Ground Control Points (GCPs) spread around a scene, we can realize this goal.The GCPs are specific pixels in an image for which the output map coordinates are known, they consist of two X, Y pairs of coordinates: i.e. • Input Coordinates: usually image coordinates in pixel unit; • Reference Coordinates: the coordinates of the map or reference image to which the source image is being registered.In this research the GCPs have been defined by using GPS navigator, shown in Table 2.

Results and Discussion
Several classic interpolation methods representing the most commonly used techniques in remotely sensing applications have been implemented in this research.Generally, the interpolation techniques fail in estimating the correct elevation values, when implemented on areas involving urban features (i.e.buildings).In our present project, the estimated errors that were produced from the implementation of the interpolation methods, specifically research's studied area was overcome, where a new procedure depending on slicing the area into two layers has been introduced.The ground-level (i.e.baselevel) has been manipulated by utilizing different interpolation techniques (i.e.IDW, Ordinary kriging, Natural neighbor), as shown in (5,6,7) figures where the DEM of the ground level by applying different interpolations techniques can be seen.
Whereas the building layer has been reclassified by assigning the GPS recorded directly to polygon points as shown in Fig. 8. Then mosaic model consider to be good techniques for merging the two layers i.e.(building and ground layer) to get the final shape of the DEM for science college of Baghdad university as shown in Fig. 9.
The geo-referenced image of the studied area can be projected on the produced DEM to represent a 3D view.

Fig. 2 :
Fig. 2: shows the image of the studied area.