This study proposes an improved algorithm that can automatically extract the optimal boundary of an arbitrary point distribution based on the convex hull model algorithm (Qian and Liu, 2007; Liu et al., 2011; Yang, 2021). In the improved algorithm, the α-shape algorithm which automatically adjusts the 2α value is used to extract the boundary in the two-dimensional mesh. Figure 4 shows the process of extracting the boundary.

Coordinate points are projected onto a two-dimensional plane along the Z-axis, maintaining the X and Y values. The neighborhood grid point detection method determines the internal and boundary grids. Hengl demonstrated that the optimal DEM grid resolution should be half the average point distance (Heng, 2006). Therefore, the size of the grid should be set to half the average distance of the points.

The vertical lines of the x and y-axes were drawn according to the boundary grid points to generate rectangular bounding boxes. The discrete points in the non-boundary grid were excluded by triangular meshing in the inner point cloud of the rectangular region, and the discrete points in the boundary grid were identified using the convex hull algorithm (He et al., 2023). Figure 5A shows an example of a minimum convex algorithm. The minimum rectangular bounding box (A-B-C-D) is first determined; subsequently, the uppermost point y max and the lowest point y min are connected. The rectangular bounding box is divided into left and right parts. On the left side of the data, the top, bottom, and leftmost points are connected to each line segment; points a and b are furthest from the outside of the line segment and connect a to y max , a to x min , b to y min , and b to x min , respectively. The same method was applied to the four lines formed by points a and b, and all boundary points of the left half were found in turn. The right half of the data is processed according to the same method. The algorithm is capable of effectively extracting the minimum convex deformation boundary; however, it cannot concave the boundary.

The proposed improved algorithm was used to extract the point cloud concave boundary. It sets 2α based on the convex boundary segment lengths. Shorter lengths indicate adjacent points, while longer lengths suggest non-adjacent points. Therefore, the quartile of each length value of the convex boundary line is counted. The value of 2α is determined as the highest quartile. Subsequently, the concave boundary points are extracted. Concave boundaries are then found where convex boundary lengths exceed 2α, using the α-shape algorithm.

As shown in Figure 5B, First, the α-shape algorithm is used when line segment P Q > 2 α . A circle with PQ ′ as its diameter is rolled counterclockwise from P to Q ( P Q ′ = 2 α , P Q ′ ⊥ P Q ). When the first point is P₁ in the 2α neighborhood, calculate the modulus of P₁Q. If P 1 Q > 2 α is still true, the next point P₂ is found according to the above method. Repeat the process until Pn is found where P n Q ≤ 2 α at which point the process stops. Finally, the optimal boundary is found by connecting all the boundary points in turn.

Spatial interpolation transforms spatially continuous points into continuous data surfaces (Mao, 2007; Tang, 2014; Huang et al., 2023). Five interpolation methods (IDW, RBF, TIN, NN, and OK) were selected for DEM construction, with a comparative analysis conducted to evaluate accuracy and factors affecting interpolation errors.

(1) Interpolation method

The IDW calculated the weighted average using the distance between the sampling point and the interpolation point as the weight. The IDW interpolation is expressed as Eq. 1 (Achilleos, 2008): f x , y = ∑ i = 1 n z i d i p ∑ i = 1 n 1 d i p , (1)

where f x , y is the elevation value of the point to be measured, n is the number of sample points, z i is the elevation value of the i t h sample point, d i is the Euclidean distance from the measured to the sample point, and p is a power with a value in the range of 1–3.

The RBF is more appropriate for computing homologous spatial interpolations. RBF established in Eq. 2: f x , y = ∑ i = 1 n α i φ x − x i 2 , (2)

where φ x − x i 2 is the radial basis function, x − x i 2 is the Euclidean distance from the point to be interpolated to the sampling point, α i is the interpolation coefficient, and n neighborhood points are searched by reference to a kd-tree (Zhou et al., 2008).

The TIN builds triangles using a sequence of points (Cao et al., 2014). However, it has the disadvantage that there may be sudden changes in the edge gradient (Li and Heap, 2014).

The NN is based on the mesh division within the sample region, and attribute values of the mesh vertices and interior are obtained using the linear interpolation method (Mao, 2007).

The OK is a geostatistical method that uses the variance function as weights to unbiased optimal estimation. The OK interpolation method uses the variation function to express the spatial variation (Anderson et al., 2006; Viswanathan et al., 2015; Hekmatnejad et al., 2019), according to Eq. 3: f x , y = ∑ i = 1 n λ i z i , (3)

where λ i is the weight coefficient. The weight coefficients can be obtained by the theoretical variogram model. The variogram model and its related parameters (nugget value, variation range, and partial base value) are key factors that affect the interpolation accuracy (Zhang, 2014). Variogram function models include spherical, Gaussian, and exponential model (Zhang et al., 2022b). It is crucial to select a reasonable theoretical variogram function model and parameters for simulation accuracy. First, semi-variance values for different lag distances are calculated using the sample data, leading to the construction of an experimental variogram. The chart illustrates the spatial correlation of the sample data as distance varies. Second, the best model describing the spatial structure of data is selected based on the experimental value. Third, the nugget, range, and sill of the theoretical variogram model are estimated using the nonlinear least squares fitting method. It is necessary to evaluate the degree of fit by observing the residual, calculating the fitting coefficient, or comparing graphs.

(2)Data density

Due to the presence of duplicate and null values in the original data, the study area is divided into grids of equal spacing. Duplicate data points within each grid were filtered to retain the mean value. The maximum sampling rate is defined as the ratio of the number of original samples to the total number of grids.

The error of DEM is constructed by different density point clouds. The training data is set down to the grid resolution under densities of 70%, 60%, 50%, 40%, 30%, 20%, and 10%. The amount of extracted data is shown in Table 2.

(3) Error evaluation

Some test data are used to verify DEM accuracy. Grid elevation values from the DEM are extracted at verification point locations. Based on the extracted grid elevation values and the corresponding verification point’s elevation values, the error index is calculated, and the error of the interpolation algorithm is evaluated.

The root mean square error (RMSE), calculated as Eq. 4, and goodness of fit (R²), calculated as Eq. 5, methods are selected as error evaluation indices. When the RMSE value is smaller or R² is closer to 1, the error of the DEM simulation is optimal. RMSE = 1 n ∑ i = 1 n Z i − Z i ^ 2 , (4) R 2 = 1 − ∑ i = 1 n Z i − Z i ^ 2 ∑ i = 1 n Z i − Z i ¯ 2 , (5) where Z i ^ is the predicted value, Z i is the truth value, n is the number of samples used to validate the test sets, and Z i ¯ is the average elevation.

In order to quantitatively analyze the roughness of the bench floor, three evaluation indexes are defined as regional area, surface slope, and roughness.

(1) Regional area

The set boundary coordinates are [ x i , y i , i = 1 , 2 ⋯ , n ], and n is the number of boundary points. The regional area is expressed as: S 1 = 1 2 ∑ i = 1 n − 1 x i + 1 + x i · y i + 1 − y i . (6)

(2) Surface slope

Assuming that the fitting slope equation could be written as Z = AX + BY + D , the coefficients A, B, and D were obtained by using the least-squares method to fit the interpolated data. The surface slope ( θ ) is calculated thus: θ = ⁡ tan arccos 1 A 2 + B 2 + 1 . (7)

(3) Surface roughness

Surface roughness refers to the degree of deviation between real and ideal surfaces in the vertical direction. The standard deviation of the vertical distance from each point to the fitting surface is calculated as S d . If S d is smaller, the roughness is smaller.

The vertical distance from each point x i , y i , z i to the fitting surface after interpolation is expressed as follows: d i = A x i + B y i + D − z i A 2 + B 2 + 1 . (8)

All points on the fitting surface d ¯ have a vertical distance average and standard deviation S d such that d ¯ = 1 K ∑ i = 1 K d i , (9) S d = ∑ i = 0 K d i − d ¯ 2 K − 1 . (10)

The overcut or undercut areas can be used for terrain indexes. These areas are determined by judging whether the distance ( d i ) from each point to the fitting surface exceeds the threshold range ∆d (∆d=0.3 m). s _grid is the size of the model grid. S c = ∑ i = 1 m d i · s g r i d   i f   d i < − ∆ d S o = ∑ i = 1 m d i · s g r i d   i f   d i ≤ ∆ d S q = ∑ i = 1 m d i · s g r i d   i f   d i > + ∆ d . (11)

By studying the variations in RMSE values at different data densities, we aim to uncover the performance of DEM interpolation methods in the Tb.175-1 area and understand the impact of data density on DEM error.

Figure 10 shows the RMSE of five interpolation methods as influenced by data density. Increasing data density leads to a gradual RMSE reduction before stabilization. Notably, RBF is less affected by lower data density (<40%). However, RMSE values exceed 0.2 m for all algorithms between 10% and 40% density. IDW is notably sensitive to sampling density in terrain fluctuating areas. Between 40% and 60% density, errors decrease notably, with NN showing the largest reduction (0.27 m–0.154 m) and TIN the smallest (0.24m–0.15 m). Above 60% density, errors plateau below 0.2 m, indicating minimal impact on DEM accuracy. This shows that the interpolation methods have less influence on DEM accuracy when data density is larger. The maximum sampling rate of TB175-1 is 74% of the whole grid due to the reduction of sampling data and sampling blind area. At the maximum sampling rate, the RMSE values for OK, RBF, IDW, TIN and NN are 0.101 m, 0.13 m, 0.144 m, 0.152 m, and 0.156 m, respectively.

The DEM errors of six test areas were compared and analyzed under the maximum sampling rate. As shown in Figure 11, OK interpolation demonstrates the highest accuracy, with an RMSE below 0.22 m and R² exceeding 0.909 for the same sample data. Conversely, IDW exhibits the largest RMSE, with three test areas (Tb.218-1, Tb.218-2, and Tb.230-2) exceeding 0.32 m and R² below 0.853. Based on the algorithm’s principle, IDW estimates cannot surpass the elevation range of sampling points, leading to significant errors in terrain with limited sampling. RBF, TIN, and NN algorithms show similar RMSE values, with consistent trends across all samples.

The overall variation trend of DEM error is positively correlated with surface roughness. For example, the RMSE of the DEM errors in Tb.218-1 considering roughness are highly significant. Therefore, when the floor is very rough, the number of samples should be increased to the highest possible extent to enhance the terrain detail retention ability and reduce interpolation errors.

The above analysis shows the following. 1) The interpolation algorithms can be ranked according to RMSE values, from smallest to largest: OK, TIN, RBF, NN, and IDW. 2) When the sampling density is low, IDW shows the highest error while RBF exhibits the lowest. For sampling rates exceeding 60%, all five interpolation methods demonstrate minimal errors, showcasing good robustness. 3) OK is least impacted by terrain features when sampling density requirements are met. Therefore, the optimal OK interpolation yields the minimum accuracy error RMSE, which is close to that of the real surface.