Rectification of Bowl-Shape Deformation of Tidal Flat DEM derived from UAV Imaging

It is necessary to periodically obtain topographic maps of the geographical and environmental characteristics of tidal flats to systemically manage and monitor them. Accurate digital elevation models (DEMs) of the tidal flats are produced while using ground control points (GCPs); however, it is both complicated and difficult to conduct GPS surveys and readings of image coordinates that correspond to these because tidal flat areas are not easy to access. The position and distribution of GCPs affect DEMs, because the entire working area cannot be covered during a survey. In this study, a least-squares height-difference (LHD) DEM matching method with a polynomial model is proposed to increase the number of DEM grids while using a presecured precise DEM to rectify the distortion and bowl effect produced by unmanned aerial vehicle (UAV) images. The most appropriate result was obtained when the translation parameters were quadratic curve polynomials with an increasing number of grids and the rotation parameters were constant. The experimental results indicated that the proposed method reduced the distortion and eliminated the error caused by the bowl effect while only using a reference DEM.


Introduction
Coastal areas are vulnerable regions that are easily affected by natural and artificial environmental changes. There is growing interest in these areas in a number of associated fields, as various forms of land (sediment) have either emerged or disappeared, owing to these changes [1,2]; a unique ecological environment among these is the intertidal flat. The intertidal flat is submerged in water during high tide and its surface is fully exposed during low tide. It also possesses significant features in terms of habitat for shellfish, water purification, maintenance of diverse species, flood control, and recreational and scenic resources [3]. The Korean intertidal flat is one of the largest tidal flats in the world. It is also an internationally protected shelter for endangered migratory birds and is considered to be one of the prosperous ecosystems inhabited by various species living in clusters [4].
It is necessary to periodically obtain topographic maps that understand the geographical and environmental characteristics of tidal flats to systemically manage and monitor them. The bottom topography of an intertidal flat in a time series manner can be used, particularly to estimate the annual or seasonal variations in the sediment budget and coastal changes [5][6][7]. Very few investigations have been carried out so far despite the importance of intertidal flats, as access to such areas is restricted. An intertidal flat consists of inflow of seawater, soil with no solidification, and deposits of mud and silt. Therefore, non-intrusive and temporally regular measurements over intertidal flats are essential.

Proposed Approach of LHD Matching
A UAV DEM built without GCPs usually contains a deformation, owing to the bowl effect. The bowl effect is probably affected by location errors in the tie-point that is used by the software for self-calibration and aerial triangulation. Therefore, a non-metric camera with inaccurate interior orientation parameters (IOPs, i.e., focal length, principal point displacement, and lens distortion) is mounted on the UAV and the exterior orientation parameters (EOPs, i.e., position and rotation of the principal point) are determined by errors that are propagated by the bundle block adjustment process between overlapping images when a camera model with initially inaccurate values is used [17][18][19].
The purpose of our study was to rectify the DEM that is produced by the overlapping UAV images and photogrammetric software using LHD matching, because it is difficult to collect GCPs via field surveys in an inaccessible area, such as a tidal flat. We used only the ground coordinates of the pre-established reference DEM, in particular, to apply the LHD. The heights of the DEM were extracted by determining same-plane positions within a grid interval in both the reference and produced DEMs, as shown in Figure 1.

Proposed Approach of LHD Matching
A UAV DEM built without GCPs usually contains a deformation, owing to the bowl effect. The bowl effect is probably affected by location errors in the tie-point that is used by the software for selfcalibration and aerial triangulation. Therefore, a non-metric camera with inaccurate interior orientation parameters (IOPs, i.e., focal length, principal point displacement, and lens distortion) is mounted on the UAV and the exterior orientation parameters (EOPs, i.e., position and rotation of the principal point) are determined by errors that are propagated by the bundle block adjustment process between overlapping images when a camera model with initially inaccurate values is used [17][18][19].
The purpose of our study was to rectify the DEM that is produced by the overlapping UAV images and photogrammetric software using LHD matching, because it is difficult to collect GCPs via field surveys in an inaccessible area, such as a tidal flat. We used only the ground coordinates of the pre-established reference DEM, in particular, to apply the LHD. The heights of the DEM were extracted by determining same-plane positions within a grid interval in both the reference and produced DEMs, as shown in Figure 1.
Without the GCPs, a bowl-shaped deformation cannot be rectified, even by applying the LHD matching method while using the seven transform parameters and reference DEM (Figure 1). We proposed a new LHD matching technique that was applied by using polynomial models of the transformation parameters as functions of the number of grid spaces between the two DEMs to overcome this problem (Equation (1)).   (1)) that specifies the relationship between the reference and UAV DEMs (translations, rotation angles, and scale). The polynomial orders were selected after testing these functions. Without the GCPs, a bowl-shaped deformation cannot be rectified, even by applying the LHD matching method while using the seven transform parameters and reference DEM ( Figure 1). We proposed a new LHD matching technique that was applied by using polynomial models of the transformation parameters as functions of the number of grid spaces between the two DEMs to overcome this problem (Equation (1)). Figure 2 shows a flowchart of the steps used in the study. LHD matching is based on the three-dimensional (3-D) similarity equation (Equation (1)) that specifies the relationship between the reference and UAV DEMs (translations, rotation angles, and scale). The polynomial orders were selected after testing these functions. A seven-parameter 3-D similarity transformation was used to express the geometric relationship between the reference DEM and UAV DEM points: where , , and are the 3-D coordinates of the reference DEM; and, , , and are the 3-D coordinates of the UAV DEM.
∆ and ∆ from Equation (3) are substituted in Equation (4), which is then expressed, as follows: A seven-parameter 3-D similarity transformation was used to express the geometric relationship between the reference DEM and UAV DEM points: is an orthogonal rotation matrix that contains the rotational angles (ω, ϕ, κ) with respect to the X, Y, and Z axes, respectively.
Equation (1) can be linearized by Taylor series expansion, as shown in Equation (2), based on the theory of Rosenholm and Torlegård (1988): where the index o indicates an initial value.
Equation (2) can be rewritten in terms of the ∆X and ∆Y portions, as in Equation (3): If the Z-axis of the reference DEM is the same as that of the UAV DEM, then (4).

The equation can be linearized as in Equation
∆X and ∆Y from Equation (3) are substituted in Equation (4), which is then expressed, as follows: (5) is substituted into the Z term in Equation (2), and the transformation parameters are finally determined by Equation (6): , which is the initial height calculated while using the initial polynomial parameters. In the Z term in Equation (1), are the slopes with respect to the X and Y directions of the two grid spaces in the produced DEM. If the value of the slope is close to zero, the accuracy of the parameters is low, as shown in Equation (6).
∆ is computed by the least squares method and iteration (Z 1 to calculate the transformation parameters. The transformation parameters are then determined until the ∆ terms are nearly zero. LHD matching uses the height difference at the same-plane locations between the two DEMs. It is important that the height difference (d i ) of any point (i) should not be largely different from those of other points. Therefore, we first determined the same-plane positions between the two DEMs within a grid interval and then calculated the height differences at those points before applying the LHD matching. Only the point where d i is within 2 σ (σ is the standard deviation) of d m can participate in the LHD matching operation to limit the gross error point: where d m is the average height difference between all of the points in the search.

Experiment and Results
The test site chosen to apply the above-mentioned method was a tidal flat area in Hampyung Bay, which is located on the west coast of the Korean Peninsula (Figure 3). The site is an embayment situated in a deep inland from the west coast and it is a unique geographic feature that is formed without large rivers. The UAV images were acquired while using a DJI Phantom 4, as shown in Table 1, at 2:20 pm on August 12, 2016, at the onset of high tide from low tide. The photogrammetric process Sensors 2020, 20, 1602 6 of 14 was implemented by the PhotoScan to obtain the DEM of the site. The camera was self-calibrated during the alignment of the image without GCPs. At this stage, corresponding points were identified on the images and matched; in addition, the position of the camera for each image from onboard global navigation satellite system data and camera calibration parameters were automatically refined and stored in the EXIF metadata of the images [25]. With the dense stereo matching of the PhotoScan, a 10 cm resolution DEM was generated by the UAV images (Figure 4). Figure 5 shows the produced UAV DEM and LiDAR DEM (resolution 1 m; obtained in 2011) used as the reference DEM. The UAV and LiDAR DEMs are both digital surface models that represent the mean sea level elevations of the reflective surfaces of all features. The difference between the two DEMs due to vegetation is not high, because most of the test area is covered by tidal flats and the data includes roads, paddy fields, and wild land.     because most of the test area is covered by tidal flats and the data includes roads, paddy fields, and wild land.    The IOPs and EOPs of the camera influence the accuracy of the DEM derived from UAV photogrammetry [26]. The EOPs were calculated by the position of the GPS equipped on the DJI Phantom 4 and bundle adjustment of the PhotoScan; however, these were not accurate. GCPs must be installed and surveyed on the tidal flat area to obtain the accurate DEM by correction of the IOPs and EOPs. However, the surveying of GCPs is almost impossible, as access to the muddy tidal flat of the Southwest coast is very difficult. Therefore, in this study, the produced DEM was corrected while using a pre-secured LiDAR DEM with a new LHD DEM matching method. In the experiment, we applied the diagonal direction (Line A in Figure 5) as the flight path of the UAV DEM to determine the polynomial models of the transformation. Next, we used the new LHD matching method with Table 2 shows the estimated constant terms (six parameters) of the polynomial model in all cases. The results of scale factor are almost 1.0 in all cases; therefore, we have omitted the scale representation. Standard deviations of shift and rotation parameters that are determined by the least squares solution do not exceed 0.3 m and 0.1°, respectively, except Cases 10 and 11. In particular, the six parameters of Case 4 and Case 8 are similar; however, Case 4 has a higher accuracy than Case 8. Overall, Cases 10 and 11 are the least reliable. The degree of the polynomial model does not necessarily have to be high or overparameterized for using it with this method; this shows that the reliability of the parameter can be decreased. The IOPs and EOPs of the camera influence the accuracy of the DEM derived from UAV photogrammetry [26]. The EOPs were calculated by the position of the GPS equipped on the DJI Phantom 4 and bundle adjustment of the PhotoScan; however, these were not accurate. GCPs must be installed and surveyed on the tidal flat area to obtain the accurate DEM by correction of the IOPs and EOPs. However, the surveying of GCPs is almost impossible, as access to the muddy tidal flat of the Southwest coast is very difficult. Therefore, in this study, the produced DEM was corrected while using a pre-secured LiDAR DEM with a new LHD DEM matching method. In the experiment, we applied the diagonal direction (Line A in Figure 5) as the flight path of the UAV DEM to determine the polynomial models of the transformation. Next, we used the new LHD matching method with the following eleven cases to obtain the transformation parameters: Case 1: 7 parameters (X 0 , Y 0 , Z 0 , ω 0 , ϕ 0 , κ 0 , s) Case 2: 10 parameters (X 0 , Y 0 , Z 0 , X 1 , Y 1 , Z 1 , ω 0 ,ϕ 0 , κ 0 , s) Case 3: 13 parameters (X 0 , Y 0 , Z 0 , X 1 , Y 1 , Z 1 , ω 0 ,ϕ 0 , κ 0 , ω 1 ,ϕ 1 , κ 1 , s) Case 4: 13 parameters (X 0 , Y 0 , Z 0 , X 1 , Y 1 , Z 1 , X 2 , Y 2 , Z 2 , ω 0 ,ϕ 0 , κ 0 , s) Case 5: 16 parameters (X 0 , Y 0 , Z 0 , X 1 , Y 1 , Z 1 , X 2 , Y 2 , Z 2 , ω 0 ,ϕ 0 , κ 0 , ω 1 ,ϕ 1 , κ 1 , s) Case 6: 19 parameters (X 0 , Y 0 , Z 0 , X 1 , Y 1 , Z 1 , X 2 , Y 2 , Z 2 , ω 0 ,ϕ 0 , κ 0 , ω 1 ,ϕ 1 , κ 1 , ω 2 ,ϕ 2 , κ 2 , s) Case 7: 22 parameters (X 0 , Y 0 , Z 0 ,X 1 , Y 1 , Z 1 , X 2 , Y 2 , Z 2 , ω 0 ,ϕ 0 , κ 0 , ω 1 ,ϕ 1 , κ 1 , ω 2 ,ϕ 2 , κ 2 , ω 3 , ϕ 3 , κ 3 , s) Table 2 shows the estimated constant terms (six parameters) of the polynomial model in all cases. The results of scale factor are almost 1.0 in all cases; therefore, we have omitted the scale representation. Standard deviations of shift and rotation parameters that are determined by the least squares solution do not exceed 0.3 m and 0.1 • , respectively, except Cases 10 and 11. In particular, the six parameters of Case 4 and Case 8 are similar; however, Case 4 has a higher accuracy than Case 8. Overall, Cases 10 and 11 are the least reliable. The degree of the polynomial model does not necessarily have to be high or overparameterized for using it with this method; this shows that the reliability of the parameter can be decreased.  The produced DEM (Figure 6a) was modified by the transformation parameters through the proposed LHD matching that was based on the eleven cases using LiDAR DEM as the reference DEM ( Figure 6b). Table 3 shows a quantitative comparison of each modified DEM with LiDAR DEM. The error before DEM correction is 1.4 m with a maximum error of 7.8 m. After DEM correction, the RMSE and maximum error fell to 1.2 m and 6.8 m, respectively, for Case 1. The accuracy slightly improved from Case 1 to Case 3; the accuracy was not improved for Case 4. The error in the DEM is relatively large when compared to DEM resolution (10 cm) due to the inflow water from the tidal flat. There is a possibility that tie-point, self-calibration, and block-adjustment errors in PhotoScan due to the reflection of the underwater area during the process of stereo matching affect the bowl effect.   Figure 7 shows a visual comparison before and after correction for Cases 1, 2, 3, 10, and 11. The left side of Figure 7 shows the point clouds from the UAV images on the mesh of the LiDAR DEM in the lateral view of the 3-D surface, and the right side shows a comparison of the deviations in height between the UAV DEMs and LiDAR DEM.

Case # Parameters
On the left side of Figure 7, it can be seen that the original UAV DEMs (Figure 7a), and the corrected Cases 1 and 2 DEMs (Figure 7b,c), are almost parabolic in shape (dashed lines), and point clouds are inconsistent with the LiDAR mesh (elliptical shape), whereas the corrected DEMs in Cases 3, 10, and 11 (Figure 7d-f) closely match the LiDAR mesh. In the color maps on the right side of Figure 7, Cases 1 and 2 DEMs are lower than the LiDAR DEM in the central part and higher on the outer part, whereas the heights of the DEMS for Cases 3, 10, and 11 are closer to that of the LiDAR DEM. As the corrected DEMs for Cases 3-11 were almost similar, these comparisons were omitted. For a detailed comparison of the corrected DEMs, Figure 8 presents the height deviation in each of the cases between the LiDAR DEM and UAV DEM before and after the correction in the Line B profile, as shown in Figure 5. The analysis of accuracy was performed while using the height deviation between the two profiles of the corrected DEMs and LiDAR DEM. For a detailed comparison of the corrected DEMs, Figure 8 presents the height deviation in each of the cases between the LiDAR DEM and UAV DEM before and after the correction in the Line B profile, as shown in Figure 5. The analysis of accuracy was performed while using the height deviation between the two profiles of the corrected DEMs and LiDAR DEM. In Figure 8, the errors in Cases 1 and 2 were reduced when compared to the original UAV DEM. However, the error due to the bowl effect still existed even after the correction by applying parameters 7 and 10. Yet, it was observed that this error was removed by applying parameter 13 to Case 3. Therefore, if the transformation parameter is applied by the LHD method, which does not take the polynomial as a flight path of the DEM without the GCP, the bowl effect cannot be eliminated. The optimal parameters were selected from among the 11 cases by calculating the RMSE in the height deviation between the two DEMs after correction (Tables 3 and 4). In particular, the In Figure 8, the errors in Cases 1 and 2 were reduced when compared to the original UAV DEM. However, the error due to the bowl effect still existed even after the correction by applying parameters 7 and 10. Yet, it was observed that this error was removed by applying parameter 13 to Case 3. Therefore, if the transformation parameter is applied by the LHD method, which does not take the polynomial as a flight path of the DEM without the GCP, the bowl effect cannot be eliminated. The optimal parameters were selected from among the 11 cases by calculating the RMSE in the height deviation between the two DEMs after correction (Tables 3 and 4). In particular, the error sharply decreases from Case 3 to approach the spatial resolution of the UAV image (10 cm) as compared to that before the correction, and the difference in accuracy from Case 3 to Case 11 is less than cm, as shown in Table 4. Applying more parameters showed no significant effect. Table 2 shows that the parameters tended to become slightly inferior as higher-order polynomials were used (Case 10 and Case 11). Therefore, in this study, it was surmised that the optimal parameter to correct the UAV-made DEM while using the reference DEM and proposed LHD matching method was Case 4, which was similar to the bowl effect geometry. It is the secondary curved polynomial with constant rotation having the principal axis in the direction of line, as shown in Figure 5. The experimental results verified that the proposed technique reduced the error of the DEM that was produced by the UAV using only the reference DEM and also eliminated the error due to the bowl effect without requiring GCPs. Table 4. Height deviations between the two DEMs before and after correction on the line B profile in Figure 5.

Correction
Before After: Case #

Conclusions
It is important to obtain a DEM of a tidal flat periodically in order to monitor its environmental condition. It is difficult to acquire GCPs in the tidal flat, although UAV images for accurate and economic DEM production are appropriate. Therefore, we proposed a polynomial LHD DEM matching method by which the position and rotation distortion and bowl effect error were removed, without the requirement of GCPs, in a tidal flat DEM produced from UAV images using a photogrammetry software. The test site chosen to apply the proposed method was a tidal flat area in Hampyung Bay that was located on the west coast of the Korean Peninsula. The experiments were conducted in 11 cases with different correction models to determine the optimal LHD polynomial model. In the experimental investigations, we proposed an optimal model, in which the translation parameters were second-order polynomials and rotation parameters were constant as the number of DEM grids increased in the flight direction. The RMSE in the height of this case was nearly 15 cm. This result was close to the DEM resolution (10 cm) from the acquired UAV image. The experimental results showed that the proposed method could reduce the error in the tidal flat DEM while using only the reference DEM and eliminated the bowl effect error. This method is very effective for correction in the DEM produced by UAV images obtained in inaccessible areas where GCP survey is difficult.
The limitation of the method that is proposed in this paper is that the accuracy of the corrected DEM is affected by the accuracy of the reference DEM. As the height accuracy of LiDAR DEM used as the reference DEM is 15-20 cm, the UAV DEM can be regarded as having similar or lower accuracy, and it is necessary to use an independent reference point to evaluate the actual accuracy of the corrected UAV DEM.