Statistical Assessment of Some Interpolation Methods for Building Grid Format Digital Bathymetric Models

Abstract

As far as the knowledge of the seabed is concerned, both for safe navigation and for scientific research, 3D models, particularly digital bathymetric models (DBMs), are nowadays of fundamental importance. This work aimed to evaluate the quality of DBMs according to the interpolation methods applied to obtain grid format 3D surfaces from scattered sample points. Other complementary factors affecting DBM vertical accuracy, such as seabed morphological complexity and surveyed points sampling density, were also analyzed by using a factorial ANOVA experimental design. The experiments were performed on a multibeam dataset provided by the Italian Navy Hydrographic Institute (IIM) with an original resolution of 1 m × 1 m grid spacing, covering a surface of 0.24 km². Six different sectors comprising different seabed morphologies were investigated. Eight sampling densities were randomly extracted from every sector, each with four repetitions. Finally, four different interpolation methods were tested, including: radial basis multiquadric function (RBMF), ordinary kriging (OK), universal kriging (UK) and Gaussian Markov random fields (GMRF). The results demonstrated that both RBMF and OK produced very accurate DBM in areas characterized by low levels of seabed ruggedness at sampling densities of only 0.0128 points/m² (equivalent grid spacing of 8.84 m). In contrast, a higher density of 0.1024 points/m² (3.13 m grid spacing) was required to produce accurate DBM in areas with more complex seabed topography. On the other hand, UK and GMRF were strongly influenced by morphology and sampling density, yielding higher vertical random errors and more prone to slightly overestimate seabed depths. In addition, sampling density and morphology were the factors that most influenced the vertical accuracy of the interpolated DBM. In this sense, the highly statistically significant influence of the interaction between sampling density and morphology on the vertical accuracy of the interpolated DBM confirms the need to perform a preliminary analysis of seabed morphological complexity in order to increase, if necessary, the number of surveyed points in cases of complex morphologies.

1. Introduction

Knowledge of the submerged part of the Earth has always been one of the most important aspects faced by the scientific community. Navigation, port works, and the exploration of marine resources are just some of the areas that require in-depth information [1]. Bathymetric surveys are carried out whenever there is a need to precisely know the morphological trend of the seabed. They are, therefore, preliminary to the realization of maritime and river works, and are indispensable for verifying the continuity of water heads and dredging volumes [2,3]. They also play a fundamental role in monitoring the evolution of beaches and coastlines [4]. Furthermore, bathymetric surveys are useful for other aims, such as geomorphological studies, habitat mapping and many more [5,6,7].

Bathymetric surveys can be performed using different types of instrumentation. Among them, those based on the transmission of acoustic signals through the water stand out, such as single-beam sonar (SBS) [8], multiple-beam sonar (MBS) [9], and side-scan sonar (SSS) [10]. This type of instrument determines the depth by considering the time interval between the emission and the return of the sound pulse generated by the echo sounder and dividing it by two. Next, the distance from the keel to the bottom is measured by multiplying half the time taken by the impulse (exit and return) by the speed of sound in the water [11]. Another methodology is that based on satellite bathymetry, or by extracting the depths using multispectral images [12]. This approach is based on the principle that the different wavelengths of light penetrate the water at different depths to be reflected from the seabed and recorded by the satellite’s multispectral sensor. This type of remote sensing technique is now widely applied through using both medium-spatial-resolution (Landsat [13]) and high-spatial-resolution satellite imagery (QuickBird [14] or Worldview 3 [15]). In these cases, the vertical accuracy achieved for depth data took values, in terms of RMSE, of 1.41 m for Landsat [13], 0.518 m for the log-band ratio method and 0.35 m for the OBRA method working on QuickBird imagery [14], or 3.71 m for WorldView 3 [15]. Furthermore, there are several studies in the literature in which electronic navigational charts (ENC) are used as a source to extract bathymetric information [16].

Depth data obtained from the different technologies are used to ensure safe navigation by producing different type of maps. As data sample densities from hydrographic sensors have increased, methods of sea-floor representation have shifted from vector-based products such as selected soundings and contours, to gridded bathymetric models [17]. These gridded models present a height value in each node, resulting in a matrix of numbers (digital elevation model, DEM). Particularly, a 3D model of the seabed is called a digital depth model (DDM) [18].

Following the set of standards for hydrographic surveys provided by the International Hydrographic Organization (IHO) [17], bathymetric survey requirements can be highly variable depending on the intended use of the final product. In fact, the required bathymetric coverage, which is based on the combination of the survey pattern and the theoretical detection area of the survey instrumentation, depends on the formal classification of the survey area. For example, the so-called “Special Order Areas”, i.e., those areas where under-keel clearance is critical for safe navigation, require very demanding minimum standards in terms of feature detection (detection of cubic features > 1 m) and maximum allowable total vertical uncertainty (lower than $\sqrt{{0.25}^2+{\left(0.0075\ast depth\right)}^2}$ meters where seabed depth is entered in meters). In this sense, it is useful to count on spatial methods able to convert scattered seabed depth points into a 3D grid model. This way of manipulating spatial information to extract new information is usually performed with a Geographic Information System (GIS), which allows the user to perform geoprocessing tasks such as data interpolation. In this regard, spatial interpolation is the process, method, or mathematical function used to estimate an unknown attribute value at unmeasured/unsampled points from measurements made at surrounding sites (known values of sampled points) [19]. Note that obtaining a full bathymetric coverage is quite expensive, so interpolation methods could be an interesting way to reduce the price (i.e., reduce sample density) in relatively smooth seabed areas while maintaining the required minimum standards.

All interpolation techniques are based on the simple concept that the closest points have more similar values [20], which means that nearby points must be at similar heights, otherwise the field will be discontinuous. This is known as Tobler’s First Law of Geography [21], a phenomenon usually known as spatial dependence, spatial continuity or autocorrelation [22].

Interpolation methods can be classified in many ways including local/global, exact/ approximate, and deterministic/geostatistical methods.

Global interpolation methods use a single function or model to fit the entire dataset. This function is then used to estimate values at any point within the data range [23]. On the other hand, local interpolation methods use a different model for each point or local neighbourhood of points. These methods use the nearby data points to estimate the value of an unknown point [24,25].

Exact interpolation methods use the known data points to derive a surface that passes exactly through all the data points [26], while approximate interpolation methods construct interpolated surfaces that approximate but do not pass through known points [27]. In this sense, exact interpolation methods provide the exact values of the function at the known data points. By contrast, approximate interpolation methods are only accurate to a certain extent at sample points, depending on the complexity of the formula used and the number of data points available [28].

Another distinction is made between deterministic and stochastic methods. Deterministic methods do not use the statistical properties of the measured points. In this sense, they are usually applied when there is sufficient knowledge of the surface to be modeled, each predicted value being the result of a deterministic function [29]. On the other hand, the stochastic methods leverage the statistical properties of the measured points, incorporating the concept of randomness and uncertainty [30,31,32,33,34].

Finally, deterministic interpolation methods are used to create gridded surfaces from scattered sampling points (actually measured points) based on either the degree of similarity (e.g., inverse distance weighting (IDW)) or the degree of smoothing (e.g., radial basis functions (RBF)). On the other hand, geostatistical or stochastic interpolation methods consider random functions, trying to model the spatial dependence between points and usually include prediction error or uncertainty [35].

It is now known that the accuracy of the obtained DEM is affected by various factors. In fact, several studies in the literature analyze which factors may have the greatest influence on DEM accuracy. In most cases, these factors can be reduced to the following: (i) the sampling density and spatial sampling distribution; (ii) the interpolation method applied to build the grid DEM; (iii) the morphological complexity or surface variability of the work area; and (iv) the vertical and planimetric accuracy of sampled points [36,37,38,39,40,41,42,43,44,45].

However, and unlike the vertical precision of digital elevation models which has been extensively studied, there are few works focused on statistically analyzing which factors may have the greatest influence on the vertical accuracy of seabed depth modeled as a grid DBM.

The aim of this study was to statistically evaluate the vertical accuracy of digital bathymetric models (DBM) obtained from MBS surveys as a function of the sampling density, the interpolation method, and the morphological complexity of the surveyed seabed. The possible interactions between these three factors were also studied.

2. Materials and Methods

This section provides a description of the investigated study area and type of datasets. In addition, the experimental design devised to achieve the objectives proposed above is presented.

Before testing various interpolation methods to generate the 3D bathymetric models, the analysis of the seabed morphology was carried out. In doing so, six different sectors were chosen, trying to represent different morphologies with variable complexity (seabed variability and roughness). This was the first source of variation of our experimental design. Eight sampling densities, each one composed of four repetitions, were randomly extracted for each sector, constituting the second source of variation. Once these first two steps were accomplished, the different interpolation methods tested in this work were applied to each repetition of sampling density and morphologically differentiated seabed sector.

It should be noted that each of the repetitions of each sampling density, due to its random extraction, corresponded to a different spatial distribution of the sampled points. Figure 1 shows a flowchart containing the different methodological steps conducted in this study. These steps will be explained in detail in the following sections.

2.1. Study Area and Dataset

The experiment was conducted on a dataset provided by the Istituto Idrografico della Marina Militare (IIMM) located in an area close to the coast of the island of Giglio (Italy) (Figure 2). The area extends within the following WGS84 geographic coordinates: West = 10°55′34.11″, East = 10°56′0.72″, South = 42°21′6.18″, North = 42°21′19.58″, covering a surface of 0.24 km². Depth values in the selected area ranged between −5.45 m and −108 m, presenting a nominal vertical accuracy of approximately 1 cm.

In particular, the dataset consisted of 240,000 depth points obtained with a multibeam sonar during a bathymetric survey conducted in 2012. The original grid spacing of the DBM was 1 m × 1 m, which according to the International Hydrographic Organization, would correspond to a Special Order 1 area [17].

This study area was chosen because it is characterized by a high level of variability of the seabed. The entire dataset was subdivided into 24 squares (100 m × 100 m) (Figure 3), each with 10,000 depth points. The decision to divide the entire dataset into 100 m × 100 m square sectors was made to ensure distinct seafloor topographies to study the effect of morphological factors on interpolation accuracy in constructing grid DBMs. In particular, of these 10,000 points per hectare, only 95% were used to extract various subsets as seeds to generate our seabed models with the different interpolation methods. The remaining 5% were set aside and used as checkpoints for subsequent validation. This means that up to 500 checkpoints were used to assess the DBM vertical accuracy of each sector, which is enough to obtain very high reliability in the estimation of systematic and random errors [37].

Among the 24 squares generated at the beginning, only six (Q2, Q3, Q4, Q8, Q15, Q23) were selected as representative of the different investigated seabed morphologies. These differentiated seabed morphologies were detected using the Morphological Variation Index (MVI) described in Section 2.2.

Starting from the sample of 9500 points located in each of the six selected sectors, eight sampling densities were randomly extracted, ranging from a minimum of 16 points (point density of 0.0016 points/m², equivalent to a grid spacing of 25 m) to a maximum of 2048 points (point density of 0.2048 points/m², grid spacing of about 2.2 m) (see

Table 1. Tested sampling densities with their corresponding equivalent grid spacing.

Sampling Density	Points	Points/m2	Equivalent Grid Spacing (m)
SD1	16	0.0016	25.00
SD2	32	0.0032	17.68
SD3	64	0.0064	12.50
SD4	128	0.0128	8.84
SD5	256	0.0256	6.25
SD6	512	0.0512	4.42
SD7	1024	0.1024	3.13
SD8	2048	0.2048	2.21

). Four randomly extracted sample-point distributions (repetitions) for each sampling density were considered. In order to have a better distribution of sample points, randomly stratified sampling was applied to select the same number of points in each square of 25 m × 25 m.

Figure 4 shows an example of some sampling densities extracted from the initial dataset of sector Q4 (sampling densities SD1, SD3, SD5 and SD7). They represent one of the four spatial distributions (repetitions) of points extracted for each sample density and sector. Note that the process of extracting the eight sampling densities was repeated four times by randomly varying the position of the sample points.

2.2. Seabed Morpholological Complexity

To quantitatively estimate the degree of seabed morphological complexity as a guide to select morphologically differentiated sectors in the experimental design, the MVI index introduced by Alcaras et al. [46] was used. This index is given by the following equation:

$$MVI=I_{xy}·S_{xy}$$

(1)

where the first term, $I_{xy}$, represents the variation of the direction of the slope (ascending or descending), while the second term, $S_{xy}$, refers to slope variation. Both terms in Equation (1) were calculated along both directions (x and y/latitude and longitude). The results of this index are then multiplied by 100 in order to scale up the values obtained.

MVI was applied to the entire dataset, thus computing 24 MVI values (i.e., as many as the grid squares initially generated). Subsequently, and from analyzing the values of MVI obtained, only six sectors were extracted, showing relatively different MVI values, thus representing significantly different seabed morphological complexities (

Table 2. MVI values of the representative six sectors extracted potentially corresponding to different seabed morphological complexities.

MVI	Sector
0.589	Q2
0.487	Q3
1.604	Q4
0.303	Q8
0.205	Q15
0.057	Q23

).

As can be seen from the values shown in Table 2, the selected sectors ranged from a flat seabed (0.057) to a very rough and variable seabed (1.604). Figure 5 shows the different morphological complexities of the six sectors finally selected. In particular, it is worth noting that the MVI index proved to be a good tool to estimate the morphological complexity of the investigated seabed, since sector Q4, having the highest MVI value, was also the one with the highest morphological variability (Figure 5).

2.3. Interpolation Methods

Spatial interpolation is defined as the procedure used to predict the value of attributes at unobserved points within a study area covered by existing observations [47].

Four different interpolation methods were considered and tested in this study. Radial basis multiquadric function (RBMF), ordinary kriging (OK) and universal kriging (UK) were implemented using ArcGIS software (version 10.8), while Gaussian Markov random fields (GMRF) was implemented using a MATLAB code freely available at https://github.com/3DLAB-UAL/dem-gmrf (accessed on 12 April 2023).

As regards the ordinary and universal kriging, among the various semi-variograms available in ArcGIS, the so-called stable semi-variogram was chosen (exponential kernel function), optimizing all the fitting parameters initially set up in ArcGIS. The semi-variogram measures the degree of spatial autocorrelation used to assign weights in kriging interpolations. The experimental semi-variogram plots one-half of the square of the differences between samples versus their distance (h) from one another (semi-variance) [48]. The semi-variance is represented by the following expression [13]:

$$\gamma (h)=\frac{1}{2N\left(h\right)}{\displaystyle \sum }_{i=1}^N{\left(z\right(x_i)-z(x_i+h\left)\right)}^2$$

(2)

where γ(h) is the semi-variance at lag h, N(h) is the number of pairs of sample points separated by a distance h, Z(x_i) is the measured value at location x_i, and Z(x_i + h) is the measured value at location x_i + h.

It is worth noting that the semi-variogram is a measure of how well two points in a dataset are correlated as a function of their separation distance. In this way, when two points are close to each other, they tend to be more highly correlated than when they are far apart [49]. Therefore, the semi-variogram allows kriging to estimate the spatial correlation structure of the dataset and use this information to make predictions of values at unsampled locations [50,51].

2.3.1. Radial Basis Multiquadric Function

As reported in [52], radial basis functions (RBF) are conceptually akin to fitting a rubber membrane through the measured sample points while minimizing the total curvature of the surface. They include different bases or kernels that significantly affect how the rubber membrane will conform to the sample points. A radial basis function approximation takes the following form:

$$s\left(x\right)={\displaystyle \sum }_{i\in I}{y}_i\phi (||x-i||),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}x\in R^d$$

(3)

where $\phi$:[0, $\infty$) → R is a fixed univariate function and the coefficients, ${(y_i)}_{i\in I}$ are real numbers, and $||x-i||$ represents the norm, being the Euclidean norm the most common choice. Mathematically, s(x) can be conceived as a linear combination of translations of a fixed function that is radially symmetric with respect to the given norm [53].

Among the various RBF that can be found in the bibliography, Franke [54] recommends the Multiquadric (RBMF) as the one providing the best results in terms of the statistical and visual evaluation of the modeled surface. In this case, the RBF kernel $\phi$ is given by the expression $\sqrt{d^2+c^2}$, where d equals the distance from the point to every interpolated node and c is the so-called smoothing factor.

2.3.2. Ordinary Kriging

Ordinary kriging (OK) is an interpolation method that can provide predictions of unknown values of a random function, which is considered as the BLUE (best linear unbiased estimator) in the sense of least variance [55]. It is also BUE (best unbiased estimator) if data respect the bell curve. When kriging is compared with deterministic interpolators, there are major differences. In fact, the former provides uncertainty assessment, anisotropy detection or methodology assumptions [34].

Kriging aims to predict the value of the random variable $z\left(x_0\right)$ at an unsampled point $x_0$ of a geographical region [56,57] by assuming the following model:

$$z(x_0)={\displaystyle \sum }_{i=1}^n{\lambda }_{i}z\left(x_i\right)$$

(4)

where ${\lambda }_i$ are the kriging weights.

The function $z\left(x_i\right)$ is composed of a deterministic component μ and a random function $\epsilon (x_i)$ [58].

$$z\left(x_i\right)=\mu +\epsilon (x_i)$$

(5)

The deterministic component is a constant value for each $x_i$ location in each area. As commented before, this stochastic method gives the optimal prediction under the presumption that the process is second-order stationery and normally distributed [31,59].

2.3.3. Universal Kriging

Universal kriging (UK) assumes the following model [60]:

$$z\left(x_i\right)=\mu \left(x_i\right)+\epsilon (x_i)$$

(6)

where $z\left(x_i\right)$ is the variable of interest, $\mu \left(x_i\right)$ is some deterministic function, and $\epsilon (x_i)$ is a term that considers random variation [61].

UK is applied when the dependent variable does not meet the criterion of second-order stationarity necessary for OK. Second-order stationarity means that the mean and variance are the same on the entire area and that the correlation between any two observations depends only on their relative position in space [62,63]. Unlike OK, where the mean $\mu$ is assumed constant over the entire region of study, UK assumes that the mean $\mu \left(x_i\right)$ depends on the spatial location, a relationship that can be approximated by a model [64].

2.3.4. Gaussian Markov Random Fields

A Markov random field (MRF), also called the undirected graphical model, is a set of random variables (seabed depth in this case) having a Markov property described by an undirected graph [65]. Each edge-linking node in the graph represents dependency between the two nodes involved. In the case of DEM interpolation, and supposing a 4-neigborhood dependency graph for DEM grid points, two types of nodes are usually assumed: (i) estimated elevation at grid points (m) and (ii) observed elevations (z) (see Figure 6). Both affect the central cell estimate depending on prior factors (F_p) and observation factors (F_o), respectively [66].

An important special case is when the random field is Gaussian, also called Gaussian Markov random fields (GMRF). Supposing a random variable z (z being seabed depth in this case), a GMRF X(z) is defined by a mean function µ(z) = E(X(z)) and a covariance function C(z, t) = Cov(X(z), X(t)). It has the property that, for every finite collection of points {s₁, …, s_p}, x ≡ (X(s₁), …, X(s_p))^T ∼ N(µ, Σ), where Σij = C(s_i, s_j).

The reader can find a detailed explanation of the procedure and the mathematical framework applied in this work in [62]. It should be noted that the GMRF mathematical framework has the advantage of making it possible to retrieve the estimated uncertainty for each interpolated elevation point and even includes break lines (terrain discontinuities) between adjacent cells to produce high-quality topographic models.

2.4. Vertical Accuracy Evaluation

To evaluate the vertical accuracy of the models obtained with the different interpolation methods applied, a technique known as true validation was undertaken [67]. It consists of determining the difference between the calculated depth and the measured one in each independent checkpoint extracted from the initial dataset.

Finally, the statistical values of these residuals are calculated, such as mean, standard deviation, minimum, maximum and root mean square error (RMSE) [37,68].

RMSE was calculated in accordance with the formula:

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{(Z_i\left(x,y\right)-Z_i^{\prime }\left(x,y\right))}^2}{N}}$$

(7)

where N is the number of the depth points, Z_i(x, y) is the interpolated depth at the location i(x, y), and $Z_i^{\prime }\left(x,y\right)$ is the measured depth (checkpoint) at the same location i(x, y).

The estimation of random and systematic errors by using RMSE and mean statistics, respectively, assumes that the errors follow a normal distribution and there are no outliers. However, this assumption is not usually true when dealing with modelling topographic surfaces [69]. In this way, robust accuracy metrics including median and normalized median absolute deviation (NMAD) [70] were also adopted in our experiment. Median and NMAD can be considered as robust alternatives to mean and RMSE, respectively. NMAD is expressed as follow:

$$NMAD=1.4826·media{n}_{i=1,2,\dots ,t}(\left|e_i-media{n}_{j=1,2,\dots ,t}\left(e_j\right)\right|)$$

(8)

where t is the number of checkpoints, $e_i=Z_i-Z_i^{\prime }$ is the i-th residual or estimation error. The constant 1.4826 is a correction parameter, making NMAD unbiased for the normal distribution of DEM errors, disregarding the abnormality induced by outliers [71].

It should be noted that the robustness of median and NMAD is at the cost of efficiency loss for the normal distribution [72].

2.5. ANOVA Test

The analysis of variance test (ANOVA) is a widely known statistical tool used to determine the influence of independent variables or factors (aggregated at different levels or groups) on a dependent variable by comparing the variance within groups with the variance among groups, and is the initial step in the analysis of factors that could affect a given dataset [73,74].

After performing the true validation vertical accuracy assessment for the four replicates taken at each sampling density tested, a full factorial ANOVA test was carried out with a single dependent variable (error statistics such as median, RMSE or NMAD) and several independent sources of variability or factors. Factorial experiments are defined when it is possible to evaluate the effect of two or more factors on the values taken by a dependent variable. In this case, the factors analyzed were sampling density (SD), seabed morphological complexity (M), and the interpolation method (IM), including the interactions between them. Figure 7 shows the factorial ANOVA test applied in our study. Q2 to Q23 represent the different sectors selected to investigate the factor seabed morphological complexity, while SD1 to SD8 refer to the different sampling densities tested. Finally, OK, UK, RBMF and GMRF are the four interpolation methods applied, and R1 to R4 are the four replicates taken at random for each sampling density.

Tukey’s mean separation post hoc test, also known as Tukey HSD (honestly significant difference) test, was performed for those factors or interactions that were indicated as statistically significant by the ANOVA test. Tukey’s test works by calculating a critical value based on the number of groups being compared and the total sample size. This critical value represents the minimum difference between two group means that is needed to conclude that the two groups are significantly different [75]. Then, for each pair of groups being compared, Tukey’s test calculates the difference between their means and compares it to the critical value. If the difference between the two means is larger than the critical value, the two groups are considered significantly different [42,76].

3. Results and Discussion

The vertical error statistics (systematic and random errors) corresponding to the true validation accuracy assessment were calculated by means of a specific code programmed in MATLAB.

3.1. Vertical Systematic Error

In order not to overextend the presentation of the results obtained, only the systematic errors obtained using a robust bias estimator such as the median were analyzed.

Table 3. Results of ANOVA test for the median statistic.

	Sum of Squares (SS)	Degree of Freedom (DG)	Mean Sum of Squares (MS)	F	p-Value
Sampling Density (SD)	1.238929	7	0.176990	10.736	<0.01
Morphology (M)	1.296660	5	0.259332	15.731	<0.01
Interpolation Method (IM)	0.439461	3	0.146487	8.886	<0.01
SD ∙ M	6.626965	35	0.189342	11.485	<0.01
SD ∙ IM	1.241455	21	0.059117	3.586	<0.01
M ∙ IM	1.807560	15	0.120504	7.310	<0.01
SD ∙ M ∙ IM	7.168532	105	0.068272	4.141	<0.01
Error	9.495654	576	0.016486

shows that all the factors that influenced the vertical bias given by the median statistical parameter (p-value < 0.05). In particular, from analyzing the F-value, it is possible to state that morphology (M) is the factor with the highest influence (F-value = 15.731). Next, we have the significant interaction between sampling density (SD) and morphology (F-factor = 11.485), and the high impact of the influence of sampling density (F-factor = 10.736). It should be noted that the interpolation method (IM), although statistically significant, was the factor that contributed the least (F-value = 8.886) towards explaining the variability of the systematic error (i.e., vertical bias).

Figure 8 depicts the general trend of the vertical bias (median of the residuals) with respect to the factor SD. From the analysis of this figure, it is possible to conclude that the higher the density of the points surveyed to generate the seabed 3D models, the lower the vertical bias observed. In particular, there was first a decrease of positive bias (seabed depth overestimation) for the first four sampling densities, and a settlement for the last three, where a depth underestimation (slightly negative bias of about 2 cm) was found. Anyway, only the bias obtained for SD1, SD2 and SD3 (lowest sampling densities) were significantly different from zero (no bias), showing a low average positive bias that, in any case, was below 10 cm. This is likely due to the fact that UK needs high sampling densities to obtain more accurate results [77]. Note that UK is prone to overestimating the actual seabed depth, since it assumes that the mean depends on spatial location [66]. As far as the influence of the sampling density and the distribution of points is concerned, studies such as [37,38,40] confirm the dependence of the model’s accuracy on both factors, thus stating the results obtained.

The general relationship between the interpolation method applied and the computed vertical systematic error is shown in Figure 9. It can be deduced that OK and RBMF were the interpolation methods that least influenced the bias of the seabed 3D model, both showing their potential as unbiased interpolators.

By contrast, UK and GMRF demonstrated to slightly overestimate (below 5 cm on average) seabed depths. It is worth noting that UK models the spatial variability of the mean value of seabed depth by means of a deterministic trend function, which ultimately helps to over-smooth the resulting 3D model. This issue could explain the slight seabed-depth overestimation found in the tested sectors.

In the same way, GMRF interpolation method counts on a tuning parameter called the tolerance parameter or relaxation factor (σ_p) [66]. This parameter may be compared to the smoothing factor of interpolation based on radial basis functions (e.g., [78]). In fact, a low σ_p value or tolerance might not be suitable for sharp changes in elevation over short distances (high complexity morphologies). In such situations, the resulting 3D surface model would likely look over-smoothed, as it is convenient to diminish the stiffness of the interpolated surface by increasing the tolerance parameter. In this work, and for the sake of simplicity, a σ_p value of 1 m was used for all sectors, which could mean that some local roughness areas were modeled with an excessively restrictive tolerance parameter (too low), which contributed to an over-smoothing of the finally interpolated 3D-surface model.

Figure 10 depicts the specific interaction between the factors SD and IM. In other words, it shows how the intensity of the bathymetric sampling interacts with the interpolation methods tested in relation to the vertical systematic error of the obtained seabed 3D models. The results shown in this figure confirm that UK and GMRF performed worse than OK and RBMF, leading to a positive bias (seabed depth overestimation) associated with lower sampling densities. This bias is practically removed by increasing the sampling density to obtain a grid spacing less than 8.8 m (SD4). On the contrary, both OK and RBMF did not display any significant bias for the full range of sampling densities tested. In this case, it is possible to see how the bias in the case of OK and RBMF is almost always equal to zero when considering all the sampling densities tested.

Finally, Figure 11 displays the results of the vertical systematic error due to interpolation according to the significant interactions between the factors morphology and the interpolation method. In this case, it is possible to observe how all the interpolation methods tested in this work presented a significantly different behavior depending on the morphology of the seabed. In particular, it should be noted that GMRF was the interpolation method most influenced by morphology, providing both a positive and negative bias for all sectors except Q2 and Q23. UK also had a positive bias, albeit more moderate than GMRF. Once again, OK and RBMF proved to be unbiased interpolators for the full range of morphological complexities tested in this work and are qualified as the most suitable to work on any morphology even with low sampling densities. The negative bias obtained in the Q4 sector is certainly due to the previously analyzed tolerance factor. As was also found in [79,80], this parameter greatly affects the error and strongly depends on the type of morphology characterizing the sector. In the case of sector Q4, since the area is characterized by a strong variability of the seabed surface, the tolerance parameter was too low and, therefore, there was an underestimation of the interpolated surface with respect to the real one. Thus, it can be seen that the influence of the seabed/terrain morphology greatly influences the accuracy of the digital models obtained. Similar results were also found in the studies conducted by Masataka [43] and Guo et al. [42], confirming the strong influence of the slope or the terrain roughness.

3.2. Vertical Random Error

As was previously done for the vertical systematic error, a single robust estimator of random error such as NMAD was used to analyze the influence of interpolation methods on the random error of the gridded seabed surface. NMAD is an estimator of scale or variability and is insensitive to the presence of outliers. Moreover, it is not sensitive to the sample size [80].

In order to obtain a complete statistical analysis, the ANOVA test was also performed on the RMSE statistic, but in order to not overextend the presentation of the results obtained, only those obtained for the NMAD are shown. The results are reported in the Supplementary Materials, where Table S1 shows the average value of vertical RMSE ± Standard deviation (four replicates) for each sample density and morphology tested. Table S2 shows the result of the ANOVA test, while Figures S1–S6 show the error trends as a function of the parameters analyzed (SD, M, MI, and the interactions between them). The results obtained for RMSE were very similar to those obtained and shown for NMAD.

The ANOVA test results for NMAD are reported in

Table 4. ANOVA test results for vertical random error estimated from NMAD.

	Sum of Squares (SS)	Degree of Freedom (DG)	Mean Sum of Squares (MS)	F	p-Value
Sampling Density (SD)	225.8777	7	32.2682	545.538	<0.01
Morphology (M)	112.4945	5	22.4989	380.374	<0.01
Interpolation Method (IM)	51.5142	3	17.1714	290.306	<0.01
SD ∙ M	76.6551	35	2.1901	37.027	<0.01
SD ∙ IM	75.2735	21	3.5845	60.600	<0.01
M ∙ IM	12.3450	15	0.8230	13.914	<0.01
SD ∙ M ∙ IM	20.1571	105	0.1920	3.246	<0.01
Error	34.0700	576	0.0591

. As in the case of vertical bias, all the factors influenced the vertical random error given by NMAD, showing p-values clearly less than 0.05. The greatest influencing factor was SD, followed by M and IM (F-values of 545.538, 380.374 and 290.306, respectively). In this case, the factors contributed to explaining the variability of the random error much more significantly than the double and triple interactions.

Figure 12 plots SD against the vertical random error (NMAD). The analysis of this Figure suggests that increasing the number of sampling points leads to a decreasing error. In fact, a very low density of sampling points reported a high random vertical error of about 1.8 m, while the highest sampling densities, SD7 and SD8, produced a very low interpolation vertical error of about 0.1 m. In this way, it is possible to conclude that the higher the density of the points surveyed, the lower the vertical random error observed. The plot of this error against the sampling density fits well with a decreasing potential function [37]. These results are in agreement with those reported by Borga et al. [81] and Lazzaro et al. [82], confirming that the greater the sampling density, the lower the impact of the interpolation method on the vertical random error; it is simply due to the shortening of the distance between surveyed and interpolated points. As confirmed also by Liu et al. [83], the initial dataset can be reduced to a certain level without significantly decreasing the vertical accuracy of the output 3D-surface model.

Figure 13 displays the relationship between the vertical random error given by NMAD and the morphology factor. As predictable, the largest error (1.3 m) corresponded to the sector Q4, the morphologically more complex seabed sector. On the other hand, the flat seabed that characterized sector Q23 obtained the lowest interpolation error. It is important to note that the trend of the vertical error follows the previously calculated MVI values shown in Table 2. These results confirm the reliability of this index to estimate morphological complexity, also confirming the previous results obtained by Alcaras et al. [46]. It also reveals the strong dependency between the accuracy of the interpolated 3D-surface model and the morphological complexity of the surveyed seabed.

The results of the vertical interpolation error (NMAD), depending on the applied interpolation method, are depicted in Figure 14. RBMF was the significantly best-performing interpolation method, achieving a mean NMAD below 0.3 m. It was closely followed by OK, which returned a mean NMAD value of 0.35 m. In contrast, both GMRF and UK performed significantly worse, obtaining mean NMAD values of 0.75 m and 0.9 m, respectively. These results agree with those obtained by Alcaras et al. [16], confirming the good performance of OK with the so-called stable semi-variograms among the group of stochastic interpolation methods. Similarly, RBMF proved to be an outstanding and efficient method within the group of deterministic interpolators. The study conducted by Wu et al. [25] confirms the results obtained in this work, confirming that OK and RBF perform the best when building a DBM according to the values recorded of RMSE and MAE. In addition, the studies published by Bello-Pineda et al. [84] and Curtarelli et al. [85] are in line with the obtained results, highlighting the excellent performance of the OK as an interpolation method. On the contrary, UK and GMRF were prone to over-smoothing the actual seabed surface, constituting a serious problem in the case of rugged topography.

Figure 15 reports NMAD versus sampling density for each of the morphologies tested. It is worth noting that the relationship between SD and NMAD fits well to a decreasing potential function of the type $NMAD=A_M{N}^{-B}$, where B is a constant to adjust and A_M is a factor (e.g., MVI) that depends on the morphological complexity of the surveyed area. In this sense, an increase in the sampling intensity leads to a decrease in the interpolation error that is modulated as a function of the seabed roughness. Sector Q4, having the highest seabed roughness, presented the largest vertical random error for SD1 (16 sample points), while sector Q23, the least morphologically complex, obtained the smallest error for SD1 (only about 0.25 m). These results confirm the findings reported by Aguilar et al. [37] and the study conducted by Alcaras et al. [46], stating the need to increase the number of sampling points (progressive sampling) in seabed areas characterized by a rugged topography.

The results corresponding to the vertical random error, as a function of the combination of SD and IM, are depicted in Figure 16. Again, a decreasing potential function is drawn when plotting NMAD versus SD, although the exponent seems to be modified according to the interpolation method applied. In fact, the negative slope of the decreasing potential function is much higher for UK and GMRF than for OK and RBMF. Furthermore, the vertical random error made by the OK and RBMF interpolators is smaller than that made by the UK and GMRF for all sampling densities. Regarding the two best methods, it can be noted that, starting from SD3, the error becomes lower (around 0.4 m), decreasing more and more as the sampling density used increases. On the contrary, UK and GRMF need a higher density of points to reduce the error, particularly from SD5 onwards. These results confirm the need to increase the sampling density to obtain more accurate results, as shown in [86].

Finally, Figure 17 shows the results corresponding to the influence of the interaction between M and IM on the vertical random error of the interpolated 3D surface. As mentioned above, a sort of clustering can be distinguished into two categories. On the one hand, OK and RBMF, the best interpolators for all tested sample densities. On the other hand, the UK and GMRF performed the worst. As expected, the maximum interpolation error was obtained with UK in the most rugged sector Q4 (mean value of 1.79 m), while the minimum error was obtained with OK and RBMF in sector Q23 (close to zero, i.e., almost no error).

As can be seen in Figure 17, all the interpolation methods turned out to be influenced by the morphology in which they were applied. For example, OK and RBMF generated less vertical error than UK and GMRF in all morphologies tested, although this difference was practically negligible in the flattest morphology Q23. These results prove the good performance of the morphological index used to estimate the topographic complexity of the seabed, since the error trend is quite similar to the MVI values shown in Table 2. It also confirms the strong dependence between the interpolation method and the morphology of the studied area, a relationship also pointed out in several works such as [36,37,42,46].

4. Conclusions

The results obtained in this study confirm that the morphological complexity of the seabed, the sampling density and the interpolation method have an important and statistically significant impact on the vertical accuracy of the interpolated digital bathymetric models (grid format). The error of these 3D-surface models was estimated using the true validation method and two robust statistical parameters: the median for the vertical systematic error and the NMAD for the vertical random error.

From the analysis of the results provided by the ANOVA test, it can be stated that sampling density and morphology are the factors that most influence the accuracy of interpolated digital bathymetric models, relegating the interpolation method to third place. In addition, the statistically significant influence of the interaction between sampling density and morphology on the vertical accuracy of the interpolated 3D surface confirms the need to perform a preliminary analysis of the morphological complexity of the seabed in order to increase, when necessary, the sampling intensity (surveyed points) in cases of complex morphologies.

The results demonstrated that both radial basis multiquadric function (RBMF) and ordinary kriging (OK) interpolators were able to build accurate grid-format bathymetric models in areas characterized by low levels of variation of seabed morphology with low sampling densities (8 m grid spacing). In contrast, a higher sampling density (3 m grid spacing 3 m or less) was required to produce accurate bathymetric models in areas characterized by high level of variation of seabed morphology. In cases with the most rugged seabeds, RBMF resulted to be the most accurate interpolation method, while universal kriging (UK) and Gaussian Markov fandom field (GMRF) appeared to be strongly influenced by morphology and sampling density, showing that they are more likely to slightly overestimate seabed depths.

MVI proved to be a very powerful index to estimate the degree of roughness of the seabed, so it can be recommended to carry out a preliminary study of the morphological complexity of the seabed to be surveyed. Since every preliminary study has its price, this first reconnaissance survey would be based on a survey with a low sample density or previously available navigation charts, which may be economically affordable. It should be noted that the IHO Standards for hydrographic surveys for “Exclusive Order area of survey” would require a very expensive survey, with up to 200% feature search and 200% bathymetric coverage.

Future developments of this work and further studies will be focused on testing other interpolation methods and datasets representative of different seabed configurations. In addition, further tests will be carried out to validate the use of MVI to quantitatively represent the variation of seabed morphology.