2.2. Airborne LiDAR Data
Airborne LiDAR technique is one of the most efficient remote-sensing methods used in terrain mapping since the last decade, and it provides an accuracy almost equivalent to the conventional terrestrial techniques in DTM generation [
16]. This technique bases on laser distance measurements integrated with GNSS/IMU (IMU-Inertial Measurement Unit) data from the aircraft platform (see
Figure 2) [
17]. After an almost automatic data capturing and processing procedure, the three-dimensional (3D) coordinates of all measured points are obtained, and then digital terrain models and/or digital surface models are produced. Briefly, processing the airborne laser scanning data typically includes a few steps with preparation, GNSS data evaluation, system calibration, calculation of the coordinates of all laser points in a geodetic coordinate system, filtering and classification of the laser points. Accordingly, the GNSS data recorded during the flight are first checked, and then the flight track is calculated with relative positioning employing the double differences of kinematic observations. System calibration bases on evaluations of the corresponding measurements in swath overlap and chosen control points on the ground, i.e., corners of parking lots or other specific constitutions having even profiles that are measured with terrestrial surveying techniques or GNSS. Following the control and calibration processes, the three-dimensional coordinates of reflected laser spots covering the area flown over are calculated. These preliminarily calculated 3D coordinates of the laser-scanned spots, constitute the point cloud, include the elevations of vegetation, man-made objects as well as bare topography surface points.
The position vector,
, of the observed scanned spot
can be principally expressed as a function of the estimated GPS/INS (INS-Inertial Navigation System) trajectory (see
Figure 2), parameterized with the GPS/INS position
and orientation parameters (roll (
), pitch (
) and yaw (
))
, respectively, having the observed range,
and the encoding angle,
, of the laser beam as they are shown in
Figure 2. The observation equation in the
e-frame (Earth-centered, Earth-fixed (ECEF) frame, realized by ITRF) is as follows [
18]:
in the equation
symbolize a quantity, which varies by time. Here, the trajectory-dependent parameters of position
and attitude
, which is parameterized by three Euler angles
,
,
, range
and the encoder angle
are considered as measurements. The lever-arm
and the bore-sight
are assumed to be constant that are obtained from calibration process, and the rotation
comes as a function of geographic position according to the reference ellipsoid and provides transformation from local geodetic (
-frame) to global geodetic system (
-frame).
In the result of Equation (1), the 3D coordinates of the final laser point cloud in the
e-frame may be required to transform to a national datum (
n-frame). Consequently, it may be necessary to apply mapping projection to 3D geodetic coordinates (
m-frame; mapping frame where the coordinates are expressed in a grid (i.e., east-, north-) coordinate system in used geodetic datum) according to users’ requirements. In certain applications, the laser point cloud heights, which are defined in the ellipsoidal height system are required to be transformed into a physical height system referring to the geoid surface (
v-frame). After transformation of point heights to the national vertical datum (from
e/
n-frame to
v-frame), the projection coordinates using a suitable map projection, e.g., UTM, TM, are generated (
m-frame). Depending on the software, the transformation from
e-frame to
n-frame and
v-frame, may not be applicable. In such cases, the necessary datum transformation of point cloud coordinates is applied externally. On the other hand, most these laser scanner data processing software offer the output option for the projection coordinates with a map projection selection [
19].
Figure 3 illustrates a sequential relation (transformations) among the “global frame” and other intermediate frames that are required for generating point cloud coordinates from the airborne laser scanning data according to users’ requirements [
18,
19].
As a standard data format, the point cloud data are archived and distributed in “LAS” format. Depending on the purpose in mapping, the laser spots in a point cloud are classified in order to generate the digital terrain, surface, and/or canopy models. In this process, an appropriate filtering algorithm is used [
16].
In the quality of the 3D coordinates of the point cloud data, various error sources interfere. Mainly imperfections in the sensors and their assembly cause observational errors and thus decrease the accuracy of generated models [
18]. Among the observed quantities, the accuracy of the laser-scanner range (
in
Figure 2b) mainly depends on the precision of time-of-flight measurements. Accordingly, the determination accuracy of the range is reported as 2 cm in average [
21]. The GNSS positioning accuracy that depends on the precision of GNSS receivers as well as the GNSS processing strategies, and the precision of used IMU sensors, play a significant role in 3D coordinate accuracy of the laser scanning points. The coordinate accuracies of reference GNSS stations on the ground, which are used in relative positioning, the distance of these points to the study area, as well as the flight altitude affect the 3D position accuracy of derived point cloud. Further information on the error budget and accuracy assessment of DEM generation with airborne laser scanning technique can be found in [
18] and [
21]. After calculating the 3D coordinates of the laser spots as a point cloud dataset, the quality of secondary products, i.e., DSMs, DTMs, DCMs (digital canopy models), are also affected by the performance of applied filtering algorithm. Higher the reliability in outlier removal and classification in the filtering process, the higher the quality of the generated products. Süleymanoğlu and Soycan [
22] provides an assessment of various filtering algorithms used for DTM generation from the airborne LiDAR point cloud, with a case study in Bergama territory. Süleymanoğlu and Soycan [
22] tested six commonly used filtering algorithms and compared them in four test areas having different landscape and topographical patterns. In conclusion, the iterative polynomial fitting (IPF) algorithm was suggested for the areas having a smooth landscape, urbanized to a certain extent and agricultural patterns. However, rather than IPF, the evaluation threshold with an expand window (ETEW) filter algorithm was reported as the best performing filter method in steep areas with dense vegetation and infrastructure [
22].
In literature, it is possible to find many studies on the assessment of DTM accuracy depending on various factors including GNSS data processing approach, position accuracy and the number of reference stations, filtering algorithm for blunder detection of point cloud data, interpolation and classification algorithms in DTM generating, etc. [
23,
24,
25]. Although it changes depending on the conditions related to the aerial laser scanning experiment, the vertical accuracy by airborne laser scanning systems in open areas is reported as approximately 0.15 m root means square error (RMSE) [
25,
26]. However, with appropriate planning of the measurements, using precise measurement equipment and rigorous data processing strategies, the accuracies of the DTMs obtained from the LiDAR techniques can be further improved [
27].
In this study, airborne LiDAR data were collected using
Optech Pegasus HA-500 mounted on a Beechcraft B200 aircraft in a test flight. The measurements were organized and conducted by the general directorate of mapping [
13]. Parameters of the laser sensor and flight technical data are summarized in
Table 1. Airborne laser scanning measurements were carried out and completed between October 23–November 6, 2014. The flight plan was prepared using flight planner software FMS planner 4.7.3 [
13]. The flight altitude was planned as 1200 m, average point density was at least 8 points per m
2 and the flight included 32 strips with an overlap of 25%. The scan width was 35°, and the scan line spacing was 580 m.
The recorded GNSS/IMU data during the flight was processed by the general directorate of mapping using POSPac MMS (Position and Orientation System Postprocessing Package (POSPac) Mobile Mapping Suite (MMS)) software [
28]. The collected GNSS/IMU data were processed and adjusted referring the two CORS-TR (continuously operating reference stations—Turkey) network stations’, AYVL (Ayvalık, latitude 39°18′41.04″N and longitude 26°41′09.96″E) and KIKA (Kırkağaç, latitude 39°06′21.24″N and longitude 27°40′19.92″E), which are in 65 km and 45 km distances, respectively, from the test site. During the postprocesses, RINEX (receiver independent exchange format) data of CORS-TR stations in 1-second sampling interval were employed [
13]. General directorate of mapping preprocessed the airborne laser scanning data using LMS (LIDAR Mapping Suite) software by OPTECH, the producer company of the used laser scanning sensor. In the tests of the produced point cloud data, 51 ground points that separated as 26 control and 25 test points, were used. In the result of the preliminary tests, carried out using control points in the area, the obtained vertical accuracy by means of RMSE is reported as ±0.07 m [
13]. The resulting unclassified point data were delivered by general directorate of mapping with their universal transversal mercator (UTM) projection coordinates (grid zone 35S) and ellipsoidal heights (
) in ITRF datum.
2.3. UAV Photogrammetry Data
In this study, a second DTM generated from the point cloud data of the General Directorate of geographic information systems was employed. As different from the first dataset, this second point cloud data were obtained with the aerial photographs from the unmanned aerial vehicle using indirect georeferencing. The horizontal coordinates of the ground control points, which were used for indirect georeferencing the point cloud, were obtained with GNSS measurements in ITRF datum. The orthometric heights of the ground control points were obtained through the leveling measurements in TUDKA99 datum. As depending on the ground control points and the indirect georeferencing of the aerial photographs, the point cloud data were generated in physical height system referring to the geoid (v-frame). Having the orthometric heights of the point cloud spots and consequently generating a second DTM in regional vertical datum was crucial for calculating the geoid heights as the difference between the two DTMs.
Aerial photogrammetry is a conventional approach widely used for 3D and 2D mapping the land parts by evaluating the photographs taken from an aerial vehicle. Similar to the aerial LiDAR technique, the photogrammetric approach is also commonly used for DTM generation in many applications. Emerging the light–weight unmanned aerial vehicles being equipped with GNSS-IMU devices and starting their use in aerial photogrammetry has significantly increased the use of this technique especially in small budget projects and applications. The fundamental theory of photogrammetric map production bases on a single image and stereo image pairs orientation. Although the fundamental outline of the theory has remained the same, the developing technologies with sensors integration opportunities led the modernization in mathematical models and calculation algorithms [
30,
31,
32]. In the modernization manner, the “structure from motion” (SfM) algorithm resolves the 3D structure employing a series of overlapping, offset images as in the stereoscopic photogrammetry. Although this algorithm follows the same basic tenets in principle, it has fundamental differences from conventional photogrammetry by means of the solution strategy. In SfM, the camera locations, the geometry of the scene, orientation parameters are solved automatically without requiring a network of ground control points with known 3D coordinates. Instead, the listed parameters are solved simultaneously using iterative bundle adjustment having high redundancy based on a set of featured points, which are automatically extracted from the multiple overlapping images [
32,
33]. The development of the SfM algorithm dates back to the 1990s and bases on computer vision technology and automatic feature-matching algorithms [
32]. It has been popularized through a number of cloud-processing software, which employ the SfM algorithm for generating point-cloud data for various purposes [
32,
34]. The advantages provided by the SfM algorithm led improvements in the quality of terrain products derived from aerial photogrammetry with UAVs [
35].
Figure 4 describes the relations between object point (
A), image frame (
i), UAV frame (
u) and geodetic coordinate system shown as the global frame (
g). Based on the described observational geometry, generating DEM grids with UAV imagery bases on a number of operational steps that involve the cooperation of different sensors including magnetic compass, barometer, CMOS camera, dual-frequency GNSS, IMU with gyroscope and accelerometer.
Figure 5 shows the workflow of the process for generating the georeferenced sparse point-cloud and finally DEM through the photogrammetric method. Accordingly, after planning and carrying out the fieldwork procedures, a sparse 3D model from triangulating correspondences between multiple images in the scene is obtained with employing the SfM algorithm. Then the 3D positions of the points, which were matched with feature tracking, are recovered with the help of georeferencing. In the final section, an optimization-based procedure is applied for calculating and triangulating a dense mesh-grid data. Resulting 3D point-cloud is filtered for removing blundered points and finally, the dense point-cloud data are interpolated in order to obtain a surface model in mesh-grid form.
There are two main approaches for georeferencing the image data [
36]. One of them is the indirect approach in aerial photogrammetry applications using UAVs. In this method, the image data are georeferenced indirectly using ground control points (GCP) having known coordinates that are visible on images. For indirect georeferencing, the aerial vehicle is not necessarily to be equipped with position and orientation instruments and this is assumed as an advantage. Accordingly, indirect georeferencing provides accurate positioning results in required geodetic datum depending on the accuracy and reference datum of the used ground control points’ coordinates. The indirect georeferencing process is formulated in Equation (2).
where
represents the ground coordinates vector of the object point (
A) in geodetic coordinate system (
g-frame),
is coordinates vector of camera sensors at the time (
) in
-frame,
is the rotation matrix of
i-frame to
g-frame,
is scale factor,
represents the object coordinates vector in
-frame (see
Figure 4).
On the other hand, the second approach is the direct georeferencing. This approach uses onboard sensors for image georeferencing. In this approach, the exterior orientation is typically measured using a GNSS set and an INS unit mounted on the aircraft body. This method requires accurate sensors time synchronization and misalignment calibration.
Direct georeferencing for obtaining ground coordinates of the object point (
A) in the geodetic coordinate system (
g-frame), relies on the Equation (3). Using this equation in the direct approach requires a priori knowledge of a number of systematic parameters (
,
) stated in the equation [
37].
where
is coordinates vector of GNSS-IMU sensors at the time (
) of exposure,
is the rotation matrix of
u-frame to
g-frame,
is the rotation matrix from of
i-frame to
u-frame,
is the object coordinates vector in
-frame,
is the position vector of camera projection center in
u-frame (see
Figure 4). The
lever-arm and
boresight parameters in
-frame are constant.
In each step of the 3D point cloud generation process in UAV photogrammetry, different factors including input data, image quality, GNSS position accuracy, employed processing parameters affect the final DEM accuracy. Considering the experimental results on accuracy assessments of DEMs generated by UAV photogrammetry imagery, it is seen that ~5–6 cm vertical accuracy in RMSE can be obtained in bare areas with moderate topography [
39].
Following the described methodology and indirect georeferencing procedure, 3D coordinates of the point cloud were derived using UAV photogrammetry by the General Directorate of geographic information systems of Turkey. The vertical accuracy of the point cloud data were reported as ~5–6 cm by the data provider institution. The dataset was provided with their universal transversal mercator (UTM) projection coordinates (grid zone 35S) in ITRF datum and orthometric heights () in Turkey national vertical control network (TUDKA99) datum. The density of the provided point cloud was 450 points per m2.
2.4. DTM Generation from Point Cloud Data for Modeling Local Geoid Surface
In the tests, the local geoid surface was derived from the differences between two DTMs, that were generated with point cloud data by aerial laser scanning
(
Section 2.2) and UAV photogrammetric imagery (
) (
Section 2.3), respectively. The point heights, obtained from
are the geometrically meaningful ellipsoidal heights (
) in ITRF datum and refer to GRS80 reference ellipsoid, whereas the
provided us the physically meaningful orthometric heights (
) in regional height datum TUDKA99. Hence the local geoid heights (
) were generated from the differences of
and
as
in meter (see
Figure 6 for the definition of geoid height
at a topography point).
With following the described procedure, based on differences of
and
, the local geoid models in grids and pointwise formats were generated. Numeric tests carried out using the grid models aimed to evaluate the performance of the produced models along the entire area which involves various land cover types. On the other hand, the pointwise local geoid models were calculated at the manually selected points at the bare topography. Thus, the numeric tests with the local geoid models in pointwise format were expected to provide more realistic and reliable results regarding the performance of the suggested geoid modeling methodology. The distribution of the manually selected 250-geoid control points in the study area is shown in
Figure 7. As will be explained in detail below, the numeric tests carried out with the local geoid models in gridwise format also show the assessment results of the two different classification methodologies used in DTM generation. Accordingly, the two DTM couples generated using both progressive densification and morphologic filter algorithms, respectively, were used in deriving the two gridwise local geoid models. The obtained geoidal heights from each DTM couple were investigated through their basic statistics in the following. Both filter algorithms, as their theories are mentioned in the following, were implemented with separate software.
After obtaining the point cloud data from the data provider institutions, in the first step of the processes, the data were filtered and the DTMs were generated accordingly. The original point cloud by aerial LiDAR measurements was given to us in divided tiles. However, the second point cloud obtained from UAV photogrammetric imagery covered the entire area and provided as a single dataset for the whole area. Since the study area is rather large and therefore processing dense point clouds as a whole requires high-capacity computer processors and lengthy processing time, we divided the UAV point cloud data into 1 km × 1 km subareas at first for the filtering process. Thereafter, the piecewise DTMs generated as separate tiles were connected in order to constitute the DTM of the whole area. The gaps, which were at the tile boundaries, were interpolated from the DTM data in neighboring map sheets. However, since it was recognized that the interpolation process was unsuccessful for providing smooth and seamless transition along the tile boundaries of the DTMs, the computation strategy was altered and the DTMs were generated using the point cloud data as a whole, after blunder removal and denoising the partitioned point cloud data. As mentioned above, the generation of DTMs was carried out using two different filtering approaches in classification. In different stages of the point cloud data processing and then the DTM generating: CloudCompare, Envi LiDAR, Global Mapper GIS software were used [
40,
41,
42]. In calculating local geoid models both in gridwise and pointwise formats by subtracting the orthometric heights of the
from the ellipsoidal heights of the
, the Surfer software by Golden Software company was employed [
43]. We selected 250 geoid control points on bare topography with homogeneous distribution over the area for generating the pointwise version of the local geoid model. While determining the locations of these scattered points set, we paid special attention to include the characteristic points of the topography that can be assumed to be remained unchanged, and thus to be measured during both the aerial LiDAR and UAV photogrammetry campaigns (
Figure 7). In this way, we aimed to use the topographic points, which are common in both datasets and hence to eliminate the effect of possible environmental changes within a short time-lag between the airborne and UAV flights in the results.
In the filtering step of the point-clouds, the noise and blunders in the point clouds were cleaned using statistical outlier removal (SOR) filter algorithm in CloudCompare software [
44]. SOR filter has a rather simple computation algorithm that bases on stochastic analysis of all points in the point cloud and removes the points, which do not satisfy the assigned threshold as the criterion, from the point cloud [
45]. In filtering, for each point in the point cloud, the spatial distance to a number of neighboring points are calculated. The average of the calculated distances is taken into account. The distribution of the differences between the calculated distances from their average is assumed as Gaussian. Then any point having a distance difference, which does not fit the distribution is removed as a blunder. Therefore, the distance between a points-couple constituted based on the assigned number of points for defining the neighborhood boundaries (let us say
points around the point in question) cannot exceed
in meter that:
where
is average distance,
= 1,2,3 is a multiplier constant that is selected depending on the required interval of confidence,
is the standard deviation of the distances. Hence, the points which are farther than
are removed. In the process, the definition of
, the number of points in the neighborhood and
, the multiplier of standard deviation, are critical parameters that determine the number of points, deleted from the dataset. Hence, the success of the filtering process. In the study, in a result of a trial-and-error-based approach, we determined and used as
= 6 and
= 1, with concerning to reach the required accuracy from the generated DTMs according to our purpose. After removing the blunders with applying distance-based statistical outlier removal filter, we ran a SOR-like filter algorithm for further denoising of the point cloud data in order to obtain higher accuracy DTM data in result [
46]. The formulation of the denoising algorithm, which was applied after distance-based outlier detection filter, is similar to the SOR filter. As different from the first step, the denoising filter takes the height differences among the points into account instead of the distances and remove the points having the height difference exceeding the threshold value (
, where
is the average value of height differences,
is multiplier constant,
is the standard deviation of height difference from their average value). In applying the denoising filter, the number of points in the neighborhood was chosen 6 and the multiplier constant was 1 too.
After removing the blunders and denoising the data, the number of points in the point cloud decreased. As an example, in the point cloud derived from UAV photogrammetric imagery, while the number of points was 23,175,329 in the original dataset, it reduced to 18,231,092 after filtering. Thus, 21% of all points were removed with filtering. After the stochastic-based filtering, which was carried out automatically using SOR and denoising filters in CloudCompare software modules, the remained points in both point cloud datasets were also visually detected once again and a few more points, which were suspected to be outliers, were removed manually using the same software. The point cloud in view of the vertical section before (a) and after (b) the automatic and manual (c) cleaning processes are shown in
Figure 8 [
46]. This visualization gives an idea of the success of the applied filters and the need for manual intervention after automatic cleaning of the data.
Following the filtering and denoising processes of the point clouds, the DTMs from aerial LiDAR and UAV photogrammetry point clouds were generated separately in raster (*.geotiff) format. In this stage, the DTM generation process was carried out through different filter algorithms and obtained results were assessed in order to clarify the difference between the applied algorithms by means of calculated geoid heights [
40,
47]. As was mentioned above, in the piecewise generation of the DTMs, the artifacts occurred along the tile boundaries on the DTM pattern (see
Figure 9). Because of these artifacts, in a second trial, the piecewise point cloud datasets were combined after filtering and denoising processes and then the DTMs were generated using point clouds as a whole in the study area. In order to optimize the computational time and computer processor usage, the point densities of the point clouds were partly decreased [
46].
In the classification for DTM generation, first, the progressive densification filter with adaptive triangulated irregular network (ATIN) algorithm was applied using Envi LiDAR software [
22]. Following the generation of DTMs from both aerial LiDAR and UAV photogrammetry point clouds using Envi LiDAR software, the same production was repeated with morphologic filter by using Global Mapper software [
22,
40].
Basically, in the progressive densification filter with the ATIN algorithm, the study area is divided into subareas and the minimum elevation is locally determined. The points with minimum elevation are assumed to be ground surface points (bare topography surface). Then triangulated irregular network (TIN) is generated using the Delaunay triangulation method. Considering the distance of a point (but except the seed point) in a Delaunay cell to the TIN surface (
) and the widest of the three angles (
) with the TIN surface (see
Figure 10), it is decided whether to keep or remove the point (
). In order to make this decision, the distance and angle at the point in question are compared with predefined threshold values. If the compared values are smaller than the threshold values, the point is assumed as the ground point. This filtering algorithm is an iterative process and each time when a point is removed from the ground points dataset, the TIN is regenerated, and the process continues until all points are classified as ground or object detail [
48].
Beside of the progressive densification, the morphologic filter that we applied in DTM generation with Global Mapper software, is also mentioned in four filtering algorithm categories given in [
48]. Morphologic filtering algorithm bases on mathematical morphology principles and mainly uses two operations, erosion (
⊖) and dilation (
⊕). According to the sequential use of these operators, closing (erosion-dilation) and opening (dilation-erosion) operations are applied. In the progressive morphologic filtering algorithm, the points having minimum height values in a window are chosen and an approximate surface as geometric locations of these points is determined. Secondary surfaces are generated by carrying out an opening operation to the initially described surface. The two surfaces are compared considering a threshold value (
) and the points below the threshold value are categorized as ground points. Not only the height difference threshold value, but also the window size selection is also critical to achieving good results in a morphologic filter. The size of the selected window can be enlarged according to the size of the largest object in the workspace in order to detect and filter the non-ground objects effectively [
49]. In practice, along the filtering process, the window size is increased gradually. Because the morphologic filter with fixed window size is capable to detect and remove the measurements for the buildings and trees from the point cloud data, however, it is difficult to detect all non-ground objects having various sizes. Increasing the window size is a solution for increasing the effectiveness of the filter. The ground points are differentiated by comparing the points’ value with the threshold. The height difference threshold can be determined considering the slope of topography in a study area. Assuming the slope as constant, the relation between the maximum height difference (
) for the terrain (
), window size (
) at
iteration and terrain slope (
) is formulated as:
and accordingly, the height difference threshold
is:
where
is the initial height difference threshold,
is the slope,
is the cell size and
is the maximum height difference threshold [
49].
Applying the addressed classification filters, the DEMs (
and
) with 1-m grid resolution covering the entire study area were generated from aerial LiDAR and UAV photogrammetry point cloud data, respectively, using Envi LiDAR (the progressive densification filter with ATIN algorithm) and Global Mapper software (progressive morphologic filtering algorithm). In the results, Global Mapper software was found more successful in extracting and removing the artificial objects in order to obtain the DTMs. The superior accuracy of DTMs, generated with the Global Mapper Software and the morphologic filter is clear in the statistics from the local geoid model test results (please see in
Section 3).
2.5. Regional Geoid Models Used in Validations
Geoid models can be determined using one of a few methods given in the literature and available to the users either as an analytical function or as a raster form [
10]. In geodetic applications the geoid height parameters (
) obtained from a suitable geoid model are used for various purposes including the transformation of GNSS ellipsoidal heights to the physical heights that are in the regional vertical datum, calculating corrections and reductions of the terrestrial observations to a reference surface, inspecting the Earth’s interior, observing and predicting Earth-mass transfer, etc. As depending on the required accuracy of height information, an appropriate geoid model is chosen and used in computations. In many engineering projects, geodetic infrastructure establishment studies and large scale mapping applications, the required geoid model accuracy is lower than 5 cm, whereas in studies such as in geographic information system establishment, archeological surveys, forest applications, a geoid model with decimeter or sub-decimeter level accuracy may satisfy the required vertical accuracy [
50]. Parallel to the advances in GNSS techniques and their use in mapping and surveying applications, many countries serve the high-accuracy regional geoid model to the practitioners and surveyors as part of their national geodetic infrastructure in order to support and increase the efficiency of satellite-based positioning techniques in geodetic and surveying applications [
51,
52]. In Turkey where the study area is located in the west of the country (
Figure 1), the accuracy of the most recent publicly available regional geoid model (Turkey Geoid 2003– TG03) is ~ 8.0 cm. This accuracy is not sufficient for the large scale maps and spatial data production studies and most of the engineering projects [
14].
TG03 is one of the geoid models that was used in this study. It was compared with the derived geoid models from generated DTMs in gridwise and pointwise formats, respectively. The TG03 model was calculated by the general directorate of mapping and released to the civilian users as the most recent official regional geoid model of Turkey in 3 arc-minute resolution grid data. This model was calculated using the remove-compute-restore (RCR) method with the least-squares-collocation (LSC). The used data included the terrestrial gravity observations distributed over the country having ~7–8 mGal accuracy and ~30 arc-second spatial resolution in modified Potsdam gravity datum [
53,
54]. In addition, the marine gravity data derived from satellite altimetry products by ERS1, ERS2, and TOPEX/POSEIDON satellite observations at the surrounding seas were employed to complete the land gravity dataset in computations. EGM96 with harmonic expansion degree
= 360° was used as a reference geopotential model for TG03 [
55]. In order to consider the high-frequency contribution of the topographic masses in the model, the DTM data having 450-m spatial resolution grids on land and dense bathymetry data in shoreline areas were employed [
53]. Involving a number of GNSS/leveling data sparsely distributed over the country in computations, the geoid heights by TG03 model fitted to the ITRF96 coordinate datum and the regional vertical datum of Turkey (TUDKA99). The geoid height values at the scattered geoid control points (see in
Figure 7) were interpolated from the TG03 grid nodes using the inverse distance weighting method and compared with calculated geoid heights derived from the local geoid models by DTMs differences [
56]. Since the TG03 model has already been served as grid-data, the comparison results of the local geoid models in grid form with the TG03 model were easily generated at grid nodes and visualized as “geoid height difference surfaces” in
Section 3.
In certain applications where decimeter level geoid model accuracy is sufficient, global geopotential models (
) may provide the required accuracy. Global geopotential models are released as spherical harmonic expansion equations with coefficients up to certain degree/order (d/o). The accuracies of these global models were significantly improved since 2000 with data contribution of Earth gravity field satellite missions, namely CHAMP, GRACE, GOCE and recently GRACE-FO [
57]. Equation (7) provides the mathematical formulation of the geoid height value by means of spherical harmonic expansion coefficients [
9].
where (
θ,
λ, r) are co-latitude, longitude and geocentric radius of the computation point, respectively;
is the major axis radius of the reference ellipsoid;
GM is the product of the gravitational constant and the mass of the Earth;
and
are the fully normalized spherical harmonic coefficients, related to the normal potential of the reference ellipsoid from a specific degree of
and order
m,
are the fully normalized associated Legendre functions. The global geopotential models with their spherical harmonic coefficients are obtained from a number of international data centers. The International Center for Global Earth Models (ICGEM) serves an almost complete list of the geopotential models, which were released since 1966 [
15].
In numeric tests of this study, one of the high-resolution geopotential models, EIGEN6C4 (
=2194) (EIGEN=”European improved gravity model of the earth by new techniques”), was evaluated and compared with the local geoid model derived from the photogrammetric techniques [
58]. The spatial resolution of the model corresponding to its maximum expansion degree is ~9 km. In the maximum expansion degree, its reported accuracy by means of RMSE of geoid height differences in Turkey is ~12 cm [
57]. In the numeric tests of the study, EIGEN6C4 (with d/o of 2194) was employed and the geoid heights (Equation (7)) from the model were calculated both at the geoid control points (
Figure 7) as well as at the grid nodes for surface-based comparisons. In calculation with the geopotential model, ICGEM calculation service was used [
15]. The calculated geoid heights referenced to the ITRF datum (GRS80 reference ellipsoid), in order that they were comparable with the local geoid models by the photogrammetric techniques.
After calculating the geoid models, another critical process is validating their accuracy using independent data and appropriate method before releasing these models to the users. Independent high-accuracy ground datasets are required to evaluate their external accuracy. In this purpose high-accuracy GNSS/leveling data, astrogeodetic deflections of vertical values, as well as sea-level observations at tide-gauge stations (for validations along the coastline) are commonly used [
12]. However, local geoid data derived from the photogrammetric techniques have not been employed for either determination or validation of the geoid models so far. In this manner, the numeric tests, which were carried out for clarifying the usability of the DTMs derived from the photogrammetric techniques in geoid modeling and validations make an original contribution to the geoid determination studies. The following chapter discusses the obtained results from the numeric tests of local geoids by aerial LiDAR and UAV photogrammetry DTMs.