A Robust Algorithm of Multiquadric Method Based on an Improved Huber Loss Function for Interpolating Remote-Sensing-Derived Elevation Data Sets

Remote-sensing-derived elevation data sets often suffer from noise and outliers due to various reasons, such as the physical limitations of sensors, multiple reflectance, occlusions and low contrast of texture. Outliers generally have a seriously negative effect on DEM construction. Some interpolation methods like ordinary kriging (OK) are capable of smoothing noise inherent in sample points, but are sensitive to outliers. In this paper, a robust algorithm of multiquadric method (MQ) based on an Improved Huber loss function (MQ-IH) has been developed to decrease the impact of outliers on DEM construction. Theoretically, the improved Huber loss function is null for outliers, quadratic for small errors, and linear for others. Simulated data sets drawn from a mathematical surface with different error distributions were employed to analyze the robustness of MQ-IH. Results indicate that MQ-IH obtains a good balance between efficiency and robustness. Namely, the performance of MQ-IH is comparative to those of the classical MQ and MQ based on the Classical Huber loss function (MQ-CH) when sample points follow a normal distribution, and the former outperforms the latter two when sample points are subject to outliers. For example, for the Cauchy error distribution with the location parameter of 0 and scale parameter of 1, the root mean square errors (RMSEs) of MQ-CH and the classical MQ are 0.3916 and 1.4591, respectively, whereas that of MQ-IH is 0.3698. The performance of MQ-IH is further evaluated by qualitative and quantitative analysis through a real-world example of DEM construction with the stereo-images-derived elevation points. Results demonstrate that compared with the classical interpolation methods, including natural neighbor (NN), OK and ANUDEM (a program that calculates regular grid digital elevation models (DEMs) with sensible shape and drainage structure from arbitrarily large topographic data sets), and two versions of MQ, including the classical MQ and MQ-CH, MQ-IH has a better ability of maximally reducing the impact of outliers, while faithfully preserving terrain features. Theoretically, MQ-IH is not a promising interpolation method, and some side effects can be found from its simulation results. For example, the contours of MQ-IH are coarser than ANUDEM in some locations of the real-world study site, and its hill shade may not strictly agree with the real-world surface at rough terrain. Furthermore, the computing cost of MQ-IH is much bigger than that of the classical interpolation methods. However, compared with the classical interpolation methods, MQ-IH has significant potential for interpolating remote-sensing-derived elevation data sets.


Introduction
At present, remote-sensing techniques have been widely adopted to capture DEM data sources, such as Light Detection and Ranging (LiDAR), Interferometric Synthetic Aperture Radar (InSAR) and photogrammetry [1][2][3]. However, due to the physical limitations of sensors, multiple reflectance, occlusions and low contrast of texture, remote-sensing-derived elevation data sets always suffer from outliers [4][5][6][7]. In a statistical sense, a spatial outlier is an observation whose value is significantly inconsistent with those of its spatial neighbors [8,9]. The existence of outliers seriously affects subsequent spatial data analyses and often distorts the inference process. In principle, if the uncertainties of sample data are known, then the outliers can be tested based on prior knowledge. However, this is a non-trivial task. It has been demonstrated that the uncertainties of sample points are highly dependent on the landscapes of the study areas [10,11], and the attributes of the scanner and the scanner geometry, such as distance and surface orientation [12,13].
Diagnostic and robust statistics are two different branches that deal with the problem of outliers. In the context of DEM construction with diagnostic statistics, outliers are first detected and removed, and then the traditional interpolation methods are used to interpolate the remaining data [4,14,15]. Compared with diagnostic statistics, robust statistics have great potential to be resilient to outliers, which is stable with respect to small changes in the data, and even large changes in the underlying data pattern cannot cause a complete failure of the procedures [16,17]. Theoretically, a perfect robust interpolation method is capable of maximally resisting the effect of outliers, while faithfully constructing terrain surfaces. In this paper, we only focus on the development of robust interpolation methods.
Presently, many promising algorithms have been proposed for generating DEMs, such as natural neighbor (NN), ordinary kriging (OK), and Multiquadric method (MQ) [18][19][20]. A series of tests indicated that MQ can be considered as an excellent interpolation method in terms of timing, storage, accuracy, visual pleasantness of the surface, and ease of implementation [21,22]. However, the classical MQ is an exact interpolator and sensitive to outliers due to its sum of squared error loss function.
In the theory of statistics, the Huber loss function has been widely applied to robust regression, which reduces the effect of outliers on parameter estimation, while treating non-outliers in a more standard way [23]. The Huber loss function considers a convex differentiable cost function, which is quadratic for small errors and linear for others. However, the classical Huber loss function does not completely avoid the effect of outliers, since its influence function is a constant for points with very large errors (i.e., outliers). Therefore, the objectives of this paper are to: (1) develop a robust algorithm of MQ method based on an Improved Huber loss function (MQ-IH) that is completely resist the impact of outliers; (2) test the robustness of MQ-IH to outliers for surface modeling with simulated data sets subject to errors with different distributions; and (3) compare the performance of MQ-IH with those of the classical interpolation methods, including natural neighbor (NN), ordinary kriging (OK), ANUDEM and the classical MQ for DEM construction with stereo-images-derived elevation data set.
The remainder of this paper is organized as follows. Section 2 presents a short literature review for diagnostic and robust statistics. Section 3 gives the principle and procedure of MQ-IH for interpolations. In Section 4, a numerical test for surface modeling with five groups of simulated data sets is first employed to assess the robustness of MQ-IH, and then a real-world example for DEM construction with stereo-images-derived elevation data set is used to compare the performance of MQ-IH with those of the classical interpolation methods. Discussion and conclusions are provided in Sections 5 and 6, respectively.

Related Works
Many methods related to outlier detection have been explored in different disciplines, like statistics and geosciences. In statistics, the existing traditional outlier detection algorithms can be classified into the following five categories, clustering-based, distribution-based, depth-based, density-based and distance-based. For clustering-based approaches, the major data cluster distributions are found and the points that do not belong to any cluster are defined as outliers [24,25]. Distribution-based approaches employ a standard distribution (e.g., Gaussian and Laplacian) to fit the data set, so that data points with heavy deviations from this distribution are defined as outliers [9]. Depth-based methods organize the data in different layers of k-dimensional convex hulls, where the data in the outer layers are defined as outliers [26,27]. Density-based approaches define the outlier with the analysis of data distribution density, which uses a local outlier factor (LOF) to compare the density of each data with those of its neighbors [28]. Thus, the data with a higher LOF is considered an outlier. Distance-based methods use a distance metric to measure the relationship between data samples and define the data as an outlier, which is obviously far away from the other data [29,30]. Generally, the aforementioned methods for detecting outliers mainly focus on low dimensional and non-spatial data sets [31].
In geosciences, the methods for spatial outlier detection can be generally classified into two categories, graphic approaches and quantitative tests. Graphic approaches to outlier detection are mostly based on visualizing spatial data [32]. Podobnikar [33] defined four classes of visual methods for assessing the accuracy of DTMs: visualizations according to spatial analytical operations based on one or multiple datasets; visualizations according to spatial statistical analysis; non-spatial visualizations; and other visualization techniques/other algorithms. Additionally, scatterplot and Moran scatterplot are two representative visualization approaches [34,35]. However, the scatterplot approaches are ineffective when sample points have huge data volumes and contain several outliers. Thus, many researchers have focused on the development of quantitative methods for outlier detection. For example, Liu, et al. [4] presented a method to find outliers mixed in remote-sensing-derived elevation points in terms of some sensitive outlier indices and locally adaptive and robust statistical criteria. Lu, et al. [27] proposed a set of graph-based algorithms to identify spatial outliers in US housing and Census data. Nevertheless, when several outliers cluster together, they could mask the effect of one another by increasing the size of residuals for other observations. This is the well-known masking and swapping problem. Yet, this problem can be overcome by the method based on Mahalanobis distance and Minimum Covariance Determinant (MCD). For example, Chen, et al. [31] used this method to detect spatial outliers with single and multiple attributes in West Nile virus data set. Giménez, et al. [36] employed this method to detect multivariate outliers in the processing of RTK positions. Nurunnabi, et al. [37] proposed a diagnostic principle component analysis (PCA) for dealing with mobile LiDAR point clouds, where the outliers are detected and removed with the robust Mahalanobis distance. It should be noted that the performance of Mahalanobis distance-based method is not assured for surface modeling of DEMs in areas with complex terrain characteristics, since it only depends on a simple linear interpolation to determine the neighborhood relationship in the process of outlier detection.
In principle, the robust interpolation methods resist the effect of outliers on surface modeling based on robust loss functions. In various disciplines including computer vision, computer graphics, computational geometry and robotics, different robust methods have been explored for plane surface fitting. For instance, Fleishman, et al. [38] proposed a forward-search-based robust moving least squares technique for reconstructing piecewise smooth surfaces. The random sample consensus (RANSAC) algorithm, presented by Fischler and Bolles [39], has been frequently optimized and generalized for decreasing the effect of outliers on planar surface fitting [40][41][42]. Debese, et al. [43] developed a robust fitting method of polynomial based on M-estimator for multibeam echosounder data interpolation. Nurunnabi, et al. [6] developed two robust versions of PCA with robust covariance [44] and Projection Pursuit [45] for planar surface fitting. However, the aforementioned methods are simple linear interpolators, which is ineffective in simulating steep slope areas.
Therefore, in order to avoid the shortcomings of the aforementioned methods, a promising robust interpolation method is developed in this paper, which combines the advantages of the nonlinear interpolation method (i.e., MQ) for dealing with complex terrain surfaces, and an improved Huber loss function for decreasing the effect of outliers on surface modeling.

Proposed Robust Algorithm
In data analysis, the essential framework of smoothing divides the observations into two components including signal and noise. Smoothing interpolation attempts to recover the coherent signal and remove noise, simultaneously. The general smoothing model can be expressed as [46], where f(xi, yi) is the value of a smoothing function to be modeled based on the sample point z(xi, yi), which includes the sample error e(xi, yi); n is the number of sample points. Multiquadric method (MQ) is an analytical interpolator of representing irregular surfaces that involve the summation of equations of quadric surfaces located at significant topographic points. In terms of MQ, the smoothing function f(x, y) is expressed as [47], where m is the degree of the polynomial; pj(x, y) and βj are the jth basis and coefficient of the polynomial; rj is the distance from the interpolated point (x, y) to the jth data point; q(r) and α, respectively, represent a basis function and its weight. For MQ, where c is a shape parameter. Its optimal value depends on the number and the distribution of data points, on the data vector and on the precision of the computation [48,49]. In order to perform interpolations with MQ (i.e., Equation (2)), the variables α and β must be computed. From Appendix A, we can obtain that, λ is a positive real constant, which determines the relative importance of loss function and the weight regularization term. In the case of noisy data, one avoids overfitting by taking a smaller λ value. It can be concluded from Appendix A that the classical MQ is inherently sensitive to outliers due to the least squares loss function in the objective function. In statistical theory, the Huber loss function has been widely applied to robust regression, which reduces the effect of outliers on parameter estimation, while treating non-outliers in a more standard way [23]. However, the Huber loss function cannot completely resist the effect of outliers due to the influence function of points with large residuals being a constant. Hence, an improved version is developed in this paper. It is expressed as follows, where s is the standard deviation of simulation residuals; c1 and c2 represent two constants, which are typically set as c1 = 2.5 and c2 = 3. This is a reasonable selection taking into account the fact that for a Gaussian distribution in statistics, there are only a few residuals bigger than 2.5 s [14]. In order to improve the robustness and efficiency of s, we used the scale estimator of Rousseeuw and Croux [50], which is expressed as, where medi(x) is the median of x. As Rousseeuw and Croux mentioned in their paper [50], this scale estimator has 50% breakdown point and 58% Gaussian efficiency, whereas median absolute deviation only has 37% Gaussian efficiency.
In this article, we assume that the spatial data sets under investigation are observations on a spatially continuous geospatial variable. These types of geospatial variables show gradual and continuous variation over space, and can be represented as a surface. Hence, if a point has a large simulation residual during surface modeling, we consider it as an outlier. Equation (11) indicates that the improved Huber loss function is null for large errors, and therefore completely avoids the impact of outliers.
Taking Equation (11) as the loss function in the objective function of MQ-IH, we can obtain α and β as follows (Appendix B), the other expressions can be found in Equations (7)(8)(9)(10).
Comparing Equation (4) with Equation (13), we can find that Equation (4) is a special case of Equation (13) with I2 being the identity matrix, and all the elements of e1 and e3 being zeroes. Namely, for the classical MQ, the residuals of all sample points belong to the set A2.
The procedure of MQ-IH ( Figure 1) for carrying out robust interpolations is as follows: (1) Solving the linear systems of the classical MQ method (i.e., Equation (4)) to obtain α and β. The parameter λ is pre-determined by the k-fold cross-validation (e.g., k = 10) technique [51]. (2) Computing the residuals of all sample points with r = f − Qα − Pβ.
In this paper, the pre-set tolerance is set as 0.01. (6) Carrying out interpolations with Equation (2).

Experimental Verification
In our experimental verification, a numerical test and a real-world test were, respectively, employed to assess the performance of MQ-IH for surface modeling and for DEM construction. In the numerical test, we validated the robustness of MQ-IH and MQ based on the classical Huber loss function (MQ-CH) to outliers for surface modeling, and these results were compared with the classical MQ. In the real-world example, we compared the performance of MQ-IH with those of the classical interpolation methods including NN, OK and ANUDEM, and the two versions of MQ, including the classical MQ and MQ-CH for DEM construction with stereo-images-derived elevation data set.

Numerical Test
A mathematical surface ( Figure 2) was employed as the test surface to analyze the robustness of MQ-IH, so that the "true" output value can be pre-determined to avoid uncertainty caused by uncontrollable data errors. The mathematical surface has the following formulation: where ri is the model error of the ith sample point (xi, yi); n is the number of sample points. In the numerical tests, the model errors ri(i = 1, 2, …, n) were, respectively, drawn from one of the following distributions: Case 1: the normal distribution N (0, 1) with the mean of 0 and the standard deviation of 1. Case 2: the contaminated normal distribution (1 − θ) × N (0, 1) + θ × N (0, 5 2 ), where θ is the contaminating proportion, which is, respectively, taken as 10%, 20% and 30%.
Case 3: the Cauchy distribution C (0, 1) with the location parameter of 0 and scale parameter of 1.
In principle, if model errors were drawn from the normal distribution N (0, 1) in Case 1, sample data would not contain outliers. If the model errors were drawn from the contaminated normal distribution in Case 2, about 100 θ percent of them would belong to N (0, 5 2 ), and some excessively large values might be included due to the very large variance. Because the Cauchy distribution in Case 3 is of such heavy tails that its mean and variance do not exist, the errors from this distribution generally include some outliers. Therefore, if the model errors were drawn from the distributions including Case 2 and Case 3, the sample data from Equation (19) would be contaminated by outliers.
The computational domain of the test surface is normalized into [-3, 3] × [-3, 3] and is orthogonally divided into 101 × 101 grid cells. Namely, the grid resolution is 0.06. 51 × 51 points were randomly collected to construct the original surface. Root mean squared error (RMSE), maximum error (MaxE) and minimum error (MinE) were adopted as the accuracy measures. RMSE is expressed as, where zi and hi represent the ith true and estimated values; n is the number of grid cells, i.e. n = 101 × 101. It should be noted that RMSE is an absolute accuracy measure, and should be as small as possible for a promising interpolation method. Simulation results (Table 1) indicate that for the normal distribution in Case 1, the classical MQ obtains the best result in terms of RMSE and MinE, which is closely followed by MQ-CH and MQ-IH. For example, the RMSEs of the classical MQ and MQ-CH are 0.2111 and 0.2116, respectively, whereas that of MQ-IH is 0.2162. This is reasonable when taking into account the fact that the least squares method is theoretically optimal for the errors with a Gaussian distribution. Hence, we can conclude that MQ-IH has high efficiency for interpolating noisy sample points with a normal distribution.
For the contaminated normal distribution in Case 2, the accuracies of all interpolation methods decrease compared with that in Case 1. Additionally, with the increase of contaminating proportion, the results of all methods become poorer, except for MaxE and MinE of MQ-IH. This indicates that MQ-IH is more stable than the classical MQ and MQ-CH. In addition, MQ-IH is always more accurate than the other two methods, except for the MinE with the contaminating proportion of 0.1. For example, the MaxE of MQ-IH ranges from 0.7643 to 1.1015, whereas those of MQ-CH and the classical MQ, respectively, range from 1.0610 to 2.0798, and 1.9612 to 3.3844. This clearly shows that MQ-IH can effectively resist the effect of outliers on surface modeling.
For the Cauchy distribution in Case 3, the errors of the classical MQ are much bigger than those of MQ-IH and MQ-CH. MQ-IH obtains better results than MQ-CH. For example, the RMSEs of MQ-CH and the classical MQ are 0.3916 and 1.4591, respectively, whereas that of MQ-IH is 0.3698. This further demonstrates the high robustness of MQ-IH to outliers.
In conclusion, MQ-IH has not only high efficiency for sample points with a normal distribution, but also high robustness for sample points subject to outliers. Furthermore, MQ-IH is more accurate than MQ-CH in terms of quantitative analysis.

Study Site and Data Capture
The study site with an area of about 0.46 km 2 is located in Shandong province, China. The terrain is mainly characterized as non-karstic with some isolated trees (Figure 3). The Gauss-Krueger projection, in which each zone is 3° of longitude in width, was adopted to transform the National Geodetic Coordinate System 1980 into rectangular Cartesian coordinates. The elevation datum was National Height Datum 1985. The elevations of this study site range from 28 m to 231 m with the mean of 78.3 m and the standard deviation of 22.6 m. The survey was carried out on 12 November 2012. Two Canon EOS 5D-Mark II cameras were mounted on the UAV for taking photos. The mean flying height of the UAV above the ground (relative flying height) was about 350 m. The photographic scale was 1:12,000 with an image of 5616 × 3744 pixels and the ground resolution of 8 cm. The side overlapping proportion of two neighbor images was about 50%, and the forward overlapping proportion was 80%. All images were processed with an automatic air triangulation (AT) mapping software termed Map-AT, developed by Chinese Academy of Survey & Mapping. After the stereo-image processing, about 6797 points were extracted. Namely, the average point space (PS) is about 0.15 m. In total, 203 checkpoints ( Figure 3) randomly sampled using RTK GNSS with vertical root mean square error (RMSE) of 2 cm and horizontal RMSE of 1 cm were used to test the interpolation accuracy. The checkpoints were selected by three surveyors under the rule of data accuracy and personal safety. Namely, the signal of RTK GNSS must be strong, and some places that are difficult to reach are not sampled.

Interpolation Methods
In addition to the three versions of MQ, including MQ-IH, MQ-CH and the classical MQ, the classical interpolation methods, including NN, OK and ANUDEM, were also employed to construct DEMs. In principle, the three classical interpolation methods have their respective features. NN and OK are, respectively, representative of exact and smoothing interpolators, which are widely adopted for remote-sensing-derived DEM construction [18,20,52,53]. NN is based on Voronoi tessellation of a discrete set of spatial points. It is very simple to perform, as no parameter should be assigned. The performance of OK is seriously determined by the semivariogram, which must be manually fitted before performing interpolations. ANUDEM is an interpolation method specifically designed for the creation of hydrologically correct digital elevation models (DEMs) [54,55]. The robustness of ANUDEM is determined by the parameter "vertical standard error". In our paper, OK was performed in the GS+ geostatistics software [56]. GS+ provides all geostatistics components, from semivariance analysis through kriging and mapping, in a single integrated software program widely praised for its flexibility and friendly interface. ANUDEM was performed in the "Spatial Analyst Tools" module of ArcGIS 10.1 with the "topo to raster" command, where the optimal value of the parameter "vertical standard error" was determined by trial and error. NN was also performed with ArcGIS 10.1.
In the first group of tests, the aforementioned interpolation methods were, respectively, employed to construct 2 m DEMs with all sample points. In order to conduct a qualitative comparison, the shaded relief maps and contour lines with 5 m interval of all interpolation methods were formed. In the second group, we tested the effect of sample density on the accuracy of interpolation methods by collecting two groups of points with different sample densities from the raw points, and the interpolation results were measured by the 203 checkpoints.

Results
NN is an exact interpolator, which can show the raw nature of sample points. It can be seen from Figure 4a that there are a considerable number of densely closed contours representing artificial depressions and erratic mounds at the locations A-J. This indicates that many erroneous values exist in the raw sample points. Labeling these points on the satellite image of this study site, we found that they are mainly located either at the low contrast and featureless regions, or at steep slope regions of rugged hills, where the automatic image matching algorithm always fails. In addition, some points are located on trees, which can be considered as outliers in the context of DEM construction. In conclusion, the raw elevation points suffer from noise and outliers.  For OK (Figure 4b), the densely closed contours at the locations B-D and F disappear. For ANUDEM with 5 m vertical standard deviation error (ANUDEM5) (Figure 4c), the densely closed contours at the locations C, D, F and J disappear. For ANUDEM with 6 m vertical standard deviation error (ANUDEM6) (Figure 4d), almost all of the densely closed contours disappear, except for the locations E and J. In addition, the small closed contours sporadically distributed in Figure 4a disappear in the maps of OK and ANUDEM. This indicates that OK and ANUDEM are less influenced by outliers than NN. The above results can be expected, as both OK and ANUDEM are smoothing interpolation methods, which can smooth noise inherent in sample points. For MQ (Figure 4e), the densely closed contours at the locations A-D and F-H disappear. For MQ-CH (Figure 4f), all densely closed contours disappear. However, some small closed contours are sporadically distributed in the upper-left. In comparison, MQ-IH is more robust to outliers than the other methods, which is demonstrated by the fact that all the densely closed contours at the locations A-J completely disappear, and are replaced by smooth contour lines with a natural appearance (Figure 4g).
The shaded relief map of NN is fragmented with many isolated pits and mounds (Figures 5a). This indicates that NN is very sensitive to outliers. The map of OK is smoother than that of NN. However, some erratic pits and artificial mounds still exist, especially at the locations A-J (Figure 5b). With the adjustment of vertical standard deviation error, ANUDEM seems to have the ability to resist outliers (Figure 4c,d). However, observing the maps of ANUDEM5 and ANUDEM6, we can find that ANUDEM possesses the robustness at the cost of losing subtle terrain features (Figure 5c,d). MQ is also influenced by the outliers, proved by the isolated pits and mounds at the locations A and D-J (Figure 5e). The effect of outliers is obviously decreased by MQ-CH (Figure 5f). However, MQ-CH has a coarse surface at the upper-left, and at other locations, its surface is so smooth that many subtle terrain features are lost. Comparatively, MQ-IH has a better ability of decreasing the effect of outliers, while retaining the subtle terrain features (Figure 5g). For example, the erosion cliff in the left bottom is kept very well. The erroneous pits and mounds at the locations A-J disappear and the surface keeps good continuity at these locations. Hence, we can conclude that MQ-IH is more robust to outliers than the classical interpolation methods for interpolating remote-sensing-derived elevation data set. Based on the checkpoints, we tested the effect of sample density on the interpolation accuracy of MQ-IH. Namely, three groups of data sets with different point spaces (PSs) were, respectively, collected from the original sampling points to construct 2 m DEMs. The PSs obtained were those shown in Table 2, i.e., from 0.15 m to 0.3 m. Table 2. Relationship between the number of sample points and average point space.

Number of Sample Points 6797 3822 1699
Point space (m) 0.15 0.2 0.3 Simulation results (Table 3) demonstrate that regardless of interpolation methods, the estimation accuracy decreases with the increasing PS. Irrespective of accuracy measures, MQ-IH is always more accurate than other interpolation methods in terms of RMSE and MinE. However, the MaxE of MQ-IH is slightly bigger than those of NN, OK, ANUDEM5 and ANUDEM6, which may be caused by the fact that some feature points at rough terrain are overly smoothed by MQ. Table 3. Accuracy comparison between different interpolation methods in the real-world example (unit: m). In order to test the computational cost of the aforementioned interpolation methods, their computing time for interpolating 2 m DEM with 6797 points was shown in Table 4. It should be noted that the computing time of OK is difficult to count due to the manual operation for fitting the semivarogram. Thus, its time is not list. Results (Table 4) indicate that NN is the fastest method. This result can be expected, as NN need not solve any system of linear equations. The two versions of ANUDEM are much faster than the three versions of MQ. This is due to the fact that ANUDEM is based on a finite difference technique to perform surface modeling, and one solution of the system of linear equations can achieve the simulation results. Furthermore, the system of ANUDEM is very sparse, so it can be efficiently solved [57]. However, MQ is a point-to-point interpolator. Namely, we must solve one system of linear equations to compute the value of one point, which brings huge computational cost for modeling the whole surface. Results also show that the classical MQ has a higher speed than MQ-CH and MQ-IH. This is because the latter two employ an iterative process to decrease the effect of outliers on surface modeling.

PS Accuracy Measures
In a word, the real-world example indicates that MQ-IH is more accurate than the classical interpolation methods for DEM construction in terms of both quantitative and qualitative analysis. However, its computational cost is much bigger than that of the classical interpolation methods.

Discussion
Recently, Chen and Li [58] presented a least absolute deviation (LAD)-based MQ (MQ-L) for decreasing the impact of outliers on DEM construction. Namely, MQ-L takes the least sum of absolute deviation as the objective function. It includes two independent procedures, i.e., the selection of some significant knots with a space-filling design and the solution of a linear system with LAD. Results indicated that MQ-L is more accurate than the classical MQ when sample points are subject to outliers. However, there are three major shortcomings of MQ-L: (1) LAD does not completely avoid the impact of outliers due to the least absolute deviation loss function; (2) the optimal number of knots is difficult to determine (in practice, it is subjectively set as the tenth of the total number of sample points); and (3) LAD has a huge memory cost, which seriously restricts its widespread applications in practice. In order to compare the performance of MQ-IH with that of MQ-L, the data sets in the numerical test of this paper were simulated by MQ-L. Results (Table 1) show that MQ-L is more robust to outliers than the classical MQ when sample points are subject to outliers (i.e., Cases 2 and 3). This conclusion is consistent with that obtained in the paper of Chen and Li [58]. However, MQ-IH is always more accurate than MQ-L, regardless of error distributions. This further validates the high robustness of MQ-IH.
It should be noted that spatial outliers in a statistical sense are not definitely erroneous data although they are highly likely to be so [31]. Instead of being simply removed, spatial outliers should be further checked with some prior knowledge to determine whether they should be retained or eliminated. For example, in the real-world example of this paper, we can superimpose outliers (i.e., observations with residuals belonging to A0 or A4 in the improved Huber loss function) on satellite images for error verification. The outliers located at a smooth terrain are definitely erroneous data. However, if outliers were located at natural cliffs, steep gorges, mounds or pits, they certainly contain important terrain features and hence should be retained. So, in the context of surface modeling with MQ-IH, the residuals of these retained outliers should belong to A2 in the improved Huber loss function.
Like the classical interpolation methods, elevation sample points and contour lines captured by any instruments or techniques can be used as the input data when performing interpolations with MQ-IH. However, compared with the classical interpolation methods, MQ-IH is more effective to remote-sensing-derived elevation data sets, since they often suffer from noise and outliers due to various reasons, such as the physical limitations of sensors, multiple reflectance, occlusions and low contrast of texture.

Conclusions
In this paper, we have developed a robust method of multiquadric method (MQ) based on an Improved Huber loss function (MQ-IH) for carrying out interpolations when spatial points are subject to outliers. The numerical test with the simulated model errors show that MQ-IH is more robust than MQ based on the classical Huber loss function (MQ-CH) and the classical MQ for surface modeling when sample points suffer from outliers, and MQ-IH has high efficiency when sample points follow a normal distribution. For example, for the Cauchy error distribution with the location parameter of 0 and scale parameter of 1, the RMSEs of MQ-CH and the classical MQ are 0.3916 and 1.4591, respectively, whereas that of MQ-IH is 0.3698. The real-world example of DEM construction with stereo-images-derived elevation data set demonstrates that MQ-IH outperforms the classical interpolation methods including NN, OK and ANUDEM, and the two versions of MQ including the classical MQ and MQ-CH. Specifically, NN is seriously influenced by outliers, which leads to the densely closed contour lines with an unnatural appearance and the fragmented shaded relief map with many artificial pits and erroneous mounds. OK is slightly better than NN. Nevertheless, some unnatural contours still exist and its map is so smooth than many terrain features like erosion cliffs are lost. ANUDEM is capable of resisting the effect of outliers with an adjustment of vertical standard deviation, but this is at the cost of losing some subtle terrain features. The classical MQ suffers from densely closed contour lines, artificial pits and erroneous mounds in the shaded relief map. MQ-CH is slightly better than MQ, but its surface is so smooth that many subtle terrain features are lost. Comparatively, MQ-IH has a better ability of maximally decreasing the effect of outliers, while faithfully keeping terrain features. When interpolating the 2 m DEM with the 6797 sample points in terms of quantitative analysis, the RMSEs of NN, OK, ANUDEM with 5 m and 6 m vertical standard deviation errors, the classical MQ and MQ-CH are 1.91 m, 1.38 m, 2.38 m, 2.36 m, 2.80 m and 1.17 m, respectively, whereas that of MQ-IH is only 0.97 m. Yet, MQ-IH has much bigger computing cost than the classical interpolation methods. In addition, MQ-IH has coarser contours than ANUDEM in the upper location of the real-world study site, and its shaded relief map may not satisfactorily conform to the real-world surface at rough terrain.
Theoretically, the well performance of MQ-IH can be attributed to the high interpolation accuracy of MQ for surface modeling and the high robustness of the improved Huber loss function to outliers. Based on MQ, surfaces with complex terrain morphologies can be accurately simulated. With the improved Huber loss function, the contribution of a point to MQ-based surface modeling is significantly determined by its sample error. Theoretically, when the error of a sample point is more than three times as big as the standard deviation of all sample errors, the point is considered as an outlier and then have no effect on surface modeling. Thus, MQ-IH has high robustness to outliers. In conclusion, MQ-IH has great potential for interpolating remote-sensing-derived elevation data sets.
In principle, a data set can contain four types of points including regular observations, vertical outliers, good leverage points and bad leverage points [8]. It should be noted that the bad leverage points also belong to outliers. In this paper, MQ-IH can only resist the effect of vertical outliers, as the improved Huber loss function is not robust to the bad leverage points. At present, the least sum of trimmed squares (LTS) has been widely accepted as a robust objective function for resisting vertical outliers and bad leverage points because of its highly statistical efficiency and high breakdown point [59]. Hence, our further work is to improve the robustness of MQ by using the LTS as the loss function for remote-sensing-derived DEM construction in study sites with different terrain morphologies, such as flat, rolling, karstic and mountainous areas.
In terms of the matrix inversion lemma, α and β can be expressed as,  

Appendix B: MQ-IH Formulation
The objective function of MQ-IH takes the improved Huber loss function (i.e., Equation (4)) as the loss function. In order to carry out interpolations with MQ-IH, the parameters α and β must be pre-determined by optimizing the objective function of MQ-IH. Let where k is the iterative time.
Based on the matrix inversion lemma, we can obtain,  