IN BATHYMETRIC SURFACES: IDW OR KRIGING?

Em superficies batimétricas: idw ou krigagem?

Abstract:

The representation of the submerged relief is very importance in diverse areas of knowledge such as Projects to build or reassess port dimensions, installation of moles, ducts, marinas, bridges, tunnels, mineral prospecting, waterways, dredging, silting control of river and lakes, and others. The depths of the aquatic bodies, indispensable for the representation of those, are obtained through the bathymetric surveys. However, the result of a bathymetric sampling is a grid of points that, for itself, it is not capable of generating directly the Digital Model of Depth (DMD), being necessary the use of interpolators. Currently, there are more than 40 available scientific methods of interpolation, each one with its particularities and characteristics. This study has the objective to analise, comparing, the efficiency of Universal Kriging (UK) and of the Inverse Distance Weighted (IDW) in the computational representation of bathymetric surfaces, varying in a decreasing way the quantity of sample points. Through the results, we can be stated the superiority of the interpolator Universal Kriging in efficiency in creating DMD with basis in the bathymetric surveys data.

Keywords:
Bathymetric surveys; Interpolators; Kriging; Inverse squared distance weighted; universal kriging

Resumo:

A representação do relevo submerso é de essencial importância em diversas áreas do conhecimento como em projetos para construção ou reavaliação de dimensões portuárias, instalação de moles, dutos, marinas, pontes, túneis, prospecção mineral, cursos de água, dragagem, controle de sedimentos de rios e lagos e outros. As profundidades dos corpos aquáticos, indispensáveis para a representação destes, são obtidas através dos levantamentos batimétricos. No entanto, o produto resultante de uma batimetria é uma malha de pontos amostrais que, por si só, não é capaz de gerar diretamente o Modelo Digital de Profundidade (MDP), sendo necessário o uso de interpoladores. Até o momento existem mais de 40 métodos de interpolação disponíveis na literatura, cada um com suas particularidades e características. Este estudo teve como objetivo analisar, comparativamente, a eficiência da Krigagem Universal (KU) e do Inverso Ponderado da Distância (IPD) na representação computacional de superfícies batimétricas, variando de forma decrescente a quantidade de pontos amostrais. Através dos resultados pode-se constatar a superioridade do interpolador Krigagem Universal quanto à eficiência em criar MDP com base nos dados de levantamentos batimétricos.

Palavras-chaves:
Levantamentos batimétricos; Interpoladores; Krigagem; inverso do quadrado da distância; krigagem universal

1.Introduction

Brazil has an extensive coast and the largest hydrographic net of the world, with rivers that stand out in depth, width and extension. This affirmation alone would justify any study related to the submerged floor.

Since the early 19^th century, navigators have tried to better understand the seafloor. At present, this study is necessary in portuary works, both in the construction and dredging of new ports; leasing of gas pipelines and transoceanic telephone cables; exploration of oil and other mineral resources; environmental preservation; research activities; follow-up of erosion or silting-up processes, and, especially, in navigation (Iho 2005; Sánchez 2010).

For analysis, preparation and introduction procedures of these studies, the use of Digital Elevation Models (DEMs) is essential. These models consist of a computational mathematical representation of the distribution of a space phenomenon occurring within a region of the land surface (Felgueiras 1998).

DEMs allow from a simple three-dimensional visualization of the floor to more complex analyses, like volume calculations and generation of slope maps (Felgueiras 1998).

The depths of water bodies, essential in the construction of DEMs of seafloor, are obtained through bathymetric surveys. In spite of the growing technological evolution, single-beam bathymetric survey is still the most used technique in the whole world (Iho 2005). This technique is carried out on board of vessels using single-beam echosounders for depth measurements at high sampling rate and GPS (Global Positioning System) receptors for differential planimetric positioning (Ferreira 2015). The outcome is a mesh of three-dimensional points that, by itself, is not able to directly generate the imaged floor surface. To build a DEM that represents such morphology, it is necessary to employ interpolation techniques to estimate the depth value of non-sampled places (Camargo 1998; Silveira 2014).

Interpolators are mathematical functions used in the construction of continuous surfaces from a set of collected points. They are used for densification of a sample that does not cover the whole interest area. Interpolation techniques are based, more frequently, on the basic geography principle that near objects tend to be more correlated than distant ones (Ferreira 2015).

Many are the interpolation methods found in the literature, each one with their peculiarities and characteristics. They are basically divided into deterministic and probabilistic interpolators (Santos 2010). Both methods are based on the similarity of near points to create a spatially continuous surface. Deterministic models make estimations from mathematical functions. Probabilistic models, besides mathematical functions, apply statistical methods, allowing besides creating spatially continuous surfaces, estimating the uncertainty of predictions (Ferreira 2015).

Among the available interpolators, the Inverse Distance Weighted (IDW) deterministic interpolator and the kriging probabilistic interpolator stand out (Carvalho e Assad 2005; Silva et al. 2008).

According Azpurua and Dos Ramos (2010), Meng et al. (2013) and Merwade et al. (2006) showed that Inverse Distance Weighting (IDW) to produce better results than geostatistical methods. Conversely, Bello-Pineda and Hernández-Stefanoni (2007) showed that the kriging method was better than IDW for mapping the bathymetry of the Yucatan submerged platform (Curtarelli et al. 2015).

In inverse distance weighted, maximum and minimum interpolated values are within the range of the sampling points (Ferreira 2015). This method determines values for not sampled points using a weighted linear combination of a set of sampled points. The weight is a function of the inverse of the distance raised to any mathematical exponent (Landim 2000; Watson 1985). This method is easy to implement, with few decisions to take regarding the model parameters. However, this interpolator is sensitive to clusters and to the presence of outliers (Ferreira 2015).

Kriging is an interpolator that can be exact or smoothed depending on the model, associated to prediction error analysis. To apply kriging assumptions must be followed (Ferreira 2015).

The main difference between kriging and other interpolation methods is in the way that weighting is attributed to the different samples. In kriging, weights are determined from a space analysis based on semivariogram. In addition, kriging normally provides unbiased and minimal variance of estimates. A great advantage of kriging over inverse distance weighted and other deterministic methods is easy production of prediction maps, like prediction errors and probabilities, in other words, kriging supplies the precision associated to each estimate (Vieira 2000). A disadvantage is the necessity of a series of decision making on data modulation, tendencies, adjustment of semivariograms and choice of neighborhood size. Thus, prior to interpolation, kriging requires a detailed geostatistical analysis of the studied phenomenon.

In spite of the vast use of these interpolators, many divergences exist on their choice and use. Studies by Kravchenco and Bullock (2003), demonstrated that kriging performs a more precise description of the spatial structure of the studied phenomenon. However, the inverse distance weighted interpolator is simpler to apply and demands less time.

Better results for kriging, when compared with the inverse distance weighted method were also noted by Tabios and Salas (1985), Laslett et al. (1987) and Warrick et al. (1988). In contrast, Kanegae Júnior et al. (2006), Wollenhaupt et al. (1994) and Gotway and Hartford (1996) demonstrated that inverse distance weighted is more efficient than kriging. Silva et al. (2008) and Souza et al. (2010) did not find great differences while comparing these methods.

Such divergences can be directly associated with the amount of sampling points. In agreement with Burrough apudCamargo (1998) when data are abundant, most interpolation methods produce basically identical results. Conversely, when data is scattered as in topobathymetric surveys, deterministic methods have limitations in the representation of spatial variability.

Therefore, the aim of the present study was to compare the efficiency of kriging and inverse distance weighted in the computational representation of bathymetric surfaces, decreasingly varying the amount of sampling points.

The purpose of this article is to make a comparison on the efficiency of Universal Kriging (KU) and the Weighted Inverse of Distance (IPD) in the computational representation of bathymetric surfaces.

1.1Inverse Distance Weighted (IDW)

The inverse distance weighted determines the values for not sampled places using a weighted linear combination of a set of sampled points. The weight is a function of inverse distance raised to any mathematical exponent (Landim 2000; Watson 1985). As a result, as distance increases the weights decrease; the decrease gets more intense, with higher exponents. The exponent value can be chosen by minimizing root mean square deviation (RMSD), obtained from cross validation (Ferreira 2015). Inverse distance weighted is calculated by the following Equation 1 described by Ferreira (2015):

$\hat{Z}\left(l_0\right)=\frac{{\displaystyle \sum _{i=1}^{n}p\left(d_i\right)\text{x Z}\left({\text{l}}_{\text{i}}\right)}}{{\displaystyle \sum _{i=1}^{n}p\left(d_i\right)}}$ (1)

where: is the estimated value for place l ₀ ; n is the number of measured values, Z(l _i ), involved in the prediction; is $p(d_i)=\frac{1}{d_i{}^{Pot}}$

the weight attributed to observation i (inverse distance function); and Pot is the mathematical exponent.

Souza (2003) affirmed that the algorithm of the inverse distance weighted is what better represents the floor surface for generation of the digital elevation model (DEM), as it has the characteristic of softening the study surface.

Another important characteristic of this method is that it allows the handling of dimension parameters of the search neighborhood, the number of neighbors to be processed in the calculation and the exponent to be employed in the distance weighted.

According to Landim (2000), with this method, the results are variable, from highly biased to in favor of points nearest to results where weight is practically the same for all near points. According to this same author, the exponent has the following effects on the estimated results: low exponents (0-2) emphasize local anomalies; whereas high exponents (3-5) soften local anomalies. Higher or equal exponents to 10 result in even estimates.

1.2Kriging

Geostatistics is based on the theory of regionalized variables. Such theory assumes that the studied phenomenon is stationary (Vieira 2000; Santos 2010). Geostatistical inference is based on the assumption of three hypotheses of stationarity: first and second order stationarity and semivariogram. First-order stationarity, according to Babak and Deutsch (2009) is that in the mean is constant in every area studied. According to Banerjee et al. (2015) second order stationarity is a less restrictive condition and exists if the mean and variance of the stochastic process are independent of location and covariance exists and is dependent only on distance h. The intrinsic hypothesis is the most used because it is less restrictive (Chilès & Delfiner 2012, Siqueira et al. 2010; Lark 2012), this means that it only requires the existence and stationarity of the semivariogram without any restriction regarding the existence of variance Finite (Vieira 2000). For geostatistical studies, second order stationarity is required (Guimarães 2004).

However, according to Santos (2010), such hypothesis cannot satisfy certain phenomena; in such cases, a less restrictive hypothesis can be used, the intrinsic hypothesis or semivariogram stationarity.

Intrinsic hypothesis assumes that Z (l) exists and does not depend on location l , and that for every Δd, the difference variance [Z (l + Δd) - Z (l) ] exists and does not depend on location l , where Z (l) corresponds to an occurrence of the studied phenomenon at point l and Δd is the distance between the successive occurrences (Guimarães 2004; Santos 2010).

The semivariogram is the most used tool in Geostatistics because it requires that only the intrinsic hypothesis is satisfied (Guimarães 2004). The experimental semivariogram was obtained from the calculation of semivariances according to Equation 2:

$\hat{\gamma }\left(\Delta d\right)=\frac{1}{2N\left(\Delta d\right)}{{\displaystyle \sum _{i=1}^{N\left(\Delta d\right)}\left[Z\left(l_i\right)-\text{Z}\left({\text{l}}_{\text{i}}+\Delta d\right)\right]}}^2$ (2)

where N(Δd) is the number of pairs of Z(l _i ) and Z(l _i + Δd) values, separated by a distance Δd. It is expected that differences {[Z(l _i ) - Z(l _i + Δd)]} decrease as Δd decreases, in other words, it is expected that nearest spatial observations have a more similar behavior between each other than more distant ones. Thus, it is expected that increases with distance (Camargo 1998).

As it can be Analyzed in equation (2), in the construction of the semivariogram, all possible pairs of data are examined. If distance Δd between two points is null, the semivariance will also be. When distance Δd is small, the points to be compared are very similar and, very correlated, soon the semivariance value is reduced (Ferreira 2015). The semivariogram graphic representation is shown according Figure 1, where the following parameters are identified:

Range: distance within which samples present themselves spatially correlated;

Sill: semivariogram value corresponding to its reach. From this point, it is considered that there is no more space dependence between samples; and,

Nugget: ideally, γ(0) = 0. However, for most of the studied phenomena there is a discontinuity of the semivariogram for smaller distances than the least distance between samples, then, γ(0) ≠ 0. In agreement with Camargo (1998) a part of this discontinuity can be attributed to measurement errors, but it is impossible to quantify whether the largest contribution comes from measurement errors or from small-scale variation unnoticed by the sampling.

Figure 1:
Representation of semivariogram’s parameters.

When the semivariogram presents identical behavior for all directions it is isotropic; otherwise, it is anisotropic. If anisotropy is detected, it must be corrected through linear transformations; since it prevents the existence of stationarity, condition necessary for accuracy in analysis and estimates for the area in study (Vieira 2000; Santos 2011).

After obtaining the experimental semivariogram, it is possible to adjust it through theoretical models (Santos 2011). It is important that the adjusted model represents the trend of in relation to Δd. Thus, estimates obtained from kriging will be more accurate and, consequently more reliable (Camargo 1998). The adjustment of the theoretical semivariogram is a very important phase, and must not be carried out automatically, since all the necessary parameters to apply kriging depend on the adjusted semivariogram model.

In the literature it is possible to find several isotropic models; these contemplate semivariograms with and without sill. Among the models without sill, the exponent model is quoted; while, among those with sill (the most common), the exponential, spherical and Gaussian models can be mentioned (Vieira 2000; Santos 2011).

Models with sill are known in Geostatistics as transitive models. Some of these models reach the sill asymptotically. For such models, the reach is randomly defined as the corresponding distance at 95% of the sill. Models without sill, as the name suggests, do not reach the sill; in other words, they keep on increasing as the distance increases. These models are used to represent phenomena that have infinite scattering capacity (Camargo 1998; Vieira 2000).

The main difference between kriging and other interpolation methods is in the way that weights (pi) are attributed to different samples. In kriging, weights are determined from a spatial analysis based on the experimental semivariogram. In addition, kriging supplies, on average, unbiased estimates and minimum variance. Another interesting characteristic of kriging is that it allows the calculation of the estimate variance; in other words, kriging supplies accuracy associated to each prediction (Camargo 1998; Vieira 2000).

According to Santos (2011) if trend is detected in the data, it is necessary to use universal kriging. In this method, trend removal is done by an adjustment of low degree polynomials. Then, the remaining analytical procedure becomes an analysis of residues. Universal kriging was proposed by Journel and Matheron to resolve a problem presented by the French National Institute of Geographic (IGN), related to the mapping of an underwater surface of evident inclination (Landim et al. 2002).

2.Materials and methods

The data used in the present study were collected in December of 2010 in a bathymetric survey of one of the main dammings of the Sao Bartolomeu stream located at the Federal University of Viçosa (UFV) in Viçosa city (MG, Brazil). The studied area has approximately 8800 m², 150 m length and 66 m width.

At collection, a single-beam echobathymeter and a couple of RTK (Real Team Kinematic) GPS receptors were used. After collection, data were processed in the Hypack 2010 software producing a file with 1414 points containing planimetric coordinates and respective depths according Figure 2. After statistical analysis of data, as recommended by Ferreira (2010), it was concluded that the survey accuracy in question is in agreement with the quality standards stipulated by DHN (Office of Hydrography and Navigation) and with the IHO (International Hydrographic Organization).

Figure 2:
Location of sampled points for GRID1 (a), GRID2 (b) and GRID3 (c).

In order to reach the goals, the original file, containing 1414 points, here denominated GRID1, was randomly divided into two other files, GRID2 (Figure 2, center), containing 706 points, and GRID3 (Figure 2, right), containing 359 points. Kriging and Inverse distance weighted were applied to GRIDS 1, 2 and 3, aiming to compare the efficiency of both interpolators in presence of many and few sampling points.

Firstly, supervised Kriging was applied, thus, depth data were submitted to an exploratory analysis. Basically, this type of analysis is based on construction and graphic interpretation, calculations and statistical interpretation. Such analysis is a very important procedure, as it allows detecting the existence of outliers and/or trends that may affect interpolation (Guimarães, 2004; Vilela, 2004).

In this study, exploratory analysis consisted in obtaining trend graphs, mean, variance, standard deviation, variation coefficient (CV), maximum value, minimum value, asymmetry, kurtosis estimation and outliers detection tests.

Subsequently, geostatistical analysis was carried out to verify the existence and, in this case, to quantify the degree of spatial dependence of the attribute in study, from the adjustment of the theoretical models to experimental semivariograms.

Semivariograms were also built for directions: N-S (0°), E- W (90°), NE-SW (45°) and SE-NW (135º). After estimation of, the obtained spatial structures were analyzed; theoretical semivariogram models which better conformed to the experimental semivariograms were built from these structures and from knowledge of the phenomenon in study.

When presence of spatial dependence was noted between data, inferences were carried out for kriging for not sampled places from the measured points, according to equation 3 (Camargo 1998; Vieira 2000; Ferreira 2015).

$\hat{Z}\left(l_0\right)={\displaystyle \sum _{i=1}^n{p}_i\text{x Z}\left({\text{l}}_{\text{i}}\right)}$ (3)

where: is the depth value estimated for location l ₀ ; n is the number of measured values, Z(l _i ), involved in the estimate; p _i are weights associated to the measured values.

For interpolation using inverse distance weighted, exponents were tested adopting the one which presented better results. A minimum of 10 and maximum of 15 sampled nearest neighboring points were used. For kriging, only points within the range of reach of the spatial dependence obtained by each attribute were considered.

In this study, evaluation of the performance of inverse distance weighting and kriging interpolators was carried out by cross validation, considering the estimates of the root-mean-square deviation (RMSD), mean error (ME), coefficient of determination (R ² ) and simple linear regression parameters between observed and predicted values, angular (a) and linear coefficient (b).

According to Santos (2011), RMSD reduces if the model adopted for the theoretical semivariogram is well chosen. In this case RMSD tends to be the same as the square root of the kriging variance. Likewise, a mean discrepancy close to zero is expected, indicating accuracy in the estimation. R ² will be best when it is the same as the unit, the same occurs for the angular coefficient (a). However, the linear coefficient (b) will be best when null. In agreement with Morillo Barragán et al. (2002) the RMSD between predicted and observed depths, in addition to the ME are given, respectively, by the following Equations 4 and 5.

$RMS=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(Z_i^{PRED}-Z_i^{OBS}\right)}^2}}{n}}$ (4)

$ME={\frac{{\displaystyle \sum _{i=1}^n\left(Z_i^{PRED}-Z_i^{OBS}\right)}}{n}}^2$ (5)

where:Z _i ^PRED and Z _i ^OBS respectively correspond to predicted and observed depths; n corresponds to the number of observed values and predicted correspondences.

3.Results and discussion

Results of the exploratory data analysis can be verified according Table 1.

Table 1:
Estimates of descriptive statistics on Depth of São Bartolomeu stream

	GRID1	GRID2	GRID3
Mean (m)	-4.18	-4.17	-4.14
Median (m)	-4.46	-4.45	-4.46
Variance (m²)	1.12	1.15	1.44
Standard deviation (m)	1.06	1.07	1.20
CV (%)	25.32	25.66	28.98
Maximum (m)	0	0	0
Minimum (m)	-5.54	-5.53	-5.53
Asymmetry	2.09	2.08	2.01
Kurtosis	5.57	8.45	7.40

It can be noticed that data present a mean variability, considering variance and sampling standard deviation values. Such variability is confirmed by the variation coefficient measurement, based on limits proposed by Warrick and Nielsen (1980), who consider: low (CV ≤ 12%); average (12% < CV < 60%) and high variability (CV ≥ 60%). It is observed, by the variation coefficient, standard deviation and variance, a higher sampling variation in GRID3, comparatively to GRIDS 1 and 2, affecting prediction.

Moreover, it is emphasized that right asymmetry presented by mean estimates, median and by the coefficient of asymmetry, highlight the form of damming of the Sao Bartolomeu Stream. The data showed the presence of some values distant from mean is noticed. They may be possible outliers; however, zero value depths correspond to the margin of a damming, not being, thus, atypical values.

Based on the exploratory analysis, the trend graph was built, using Arcgis 10 software according Figure 3.

Figure 3
Graphics showing second order trend present in GRID1 (a), GRID2 (b) and GRID3 (c).

In all GRIDS, the presence of second order trend is noticed in depth data, seen in parabolas exposed in vertical plans, which is in agreement with authenticity, as it is a reservoir with intense inclination.

When trend presence is noticed proper geostatistical interpolation is applied. In view of this data characteristic, universal kriging (UK) is chosen. According to Santos (2011), UK applies an adjustment of low-degree polynomials for trend removal, allowing working with residues.

Aiming at verifying the existence of anisotropy, semivariograms were calculated for directions: N-S (0°), E-W (90°), SW-NE (45°) and NW-SE (135º), according Figure 4. It is worth mentioning that as universal kriging was chosen, the built semivariograms here correspond to residual semivariograms.

Figure 4:
Directional experimental semivariograms and adjusted direction models: N-S (0°), E-W (90°), SW-NE (45°) and NW-SE (135º) for GRID1.

By analyzing the semivariograms of GRID1 shown in Figure 4, it can be noticed that the depth variable presents practically identical spatial dependence standards up to range, that is, it presents the same spatial variability in all directions. Thus, it is concluded that the phenomenon is isotropic. Hence, a single semivariogram representing all directions can be used, entitled omnidirectional semivariogram.

Therefore, semivariogram adjustment was carried out using Arcgis 10 software. An omnidirectional semivariogram was obtained that represents the trend of in relation to Δd. The same analysis was carried out for GRIDS 2 and 3, where anisotropy presence was not detected.

The chosen theoretical models for each GRID are summarized in Table 2.

The theoretical model which better adjusted to an experimental GRID1semivariogram was the stable. In this model, it is necessary to define a parameter, which varies from 0 to 2, where the null value makes the stable model identical to the exponential model. If the parameter is defined as 2, the model becomes Gaussian. The parameter value of the stable model defined in this study was 1.432227.

Table 2:
Estimates of the variogram analysis

UK	Model	Nugget (m²)	Sill (m²)	Range (m)
GRID1	Stable	0.000	0.520	38.388
GRID2	Gaussian	0.017	0.390	23.570
GRID3	Spherical	0.037	0.666	53.497

The omnidirectional experimental semivariogram and the adjusted model can be seen according Figure 5 for three grids.

Figure 5:
Omnidirectional experimental semivariogram and adjusted model for GRID1 (a), GRID2 (b) and GRID3 (c).

After obtaining the semivariogram, universal kriging interpolation can be applied. Prior to interpolation, cross validation is carried out, allowing evaluating the performance of interpolators; however the outcomes will be presented throughout the text.

As previously mentioned, universal kriging was used to estimate points in not sampled locations.

For interpolation using inverse distance weighted, the number of neighbors to be used in the interpolation was set firstly. A minimum of 10 and maximum of 15 sampled nearest points was adopted. Subsequently, a study was carried out to define the value of the exponent used as weight. This value was chosen by analyzing several factors, such as the area characteristics and the RMSD value obtained in the cross validation, as suggested by Ferreira (2015). Exponents were tested varying from 1 to 5. Results are shown according Figure 6.

Figure 6:
Graphic representation of Exponent x RMSD analysis.

By analyzing Figure 6, the choice of higher exponents for both GRIDS becomes obvious; however it is necessary to be careful with such choice. As reported by Landim (2000), the exponent choice has the following effects on estimated results: low exponents point out local anomalies, whereas high exponents soften local anomalies. In other words, the exponent controls the importance of points around the estimated value, that is, higher exponents result in fewer distant point influences.

It was noticed that lower exponents, besides pointing out local anomalies, provide a smoother surface, this fact is explained by higher weight given to most distant points. The highest exponent, value 5, despite providing a lower RMSD value, around 0.290, for GRID1, 0.428 for GRID2 and 0.535 for the GRID 3, provides a more detailed surface; in other words, less soft. Such fact is due to higher emphasis given to the nearest points.

In view of that, exponent value 2 was chosen, mainly due to the studied area characteristics along with the vast use of this exponent in the literature (Morillo Barragán et al., 2002; Silva et al., 2008; Souza et al., 2010). It is worth pointing out that when exponent value 2 is chosen, inverse distance weighted is then called inverse squared-distance weighted (ISDW).

In this study the evaluation of the performance of ISDW and UK interpolators was carried out by cross validation. Results are shown according Table 3.

Table 3:
Presentation of main cross validation measures

	GRID1 - 1414 points		GRID2 - 706 points		GRID3 - 359 points
	UK	ISDW	UK	ISDW	UK	ISDW
RMSD (m)	0.123	0.369	0.194	0.507	0.367	0.675
ME (m)	-0.002	0.061	0.003	0.092	0.008	0.154
R 2	0.986	0.903	0.966	0.821	0.902	0.772
a (m)	0.995	0.911	0.976	0.884	0.985	0.784
b (m)	-0.018	-0.393	-0.108	-0.514	-0.067	-0.975

When analyzing Table 3, GRID by GRID, it can be noticed, through all adopted decision parameters, that UK favored higher accuracy, in both GRIDS, when compared to ISDW, fact justified by RMSD and ME values. In addition, according to Vieira (2000), simple linear regression between observed and predicted values must present R ² quite near the unit, as well as the regression coefficient "a" and intercept "b" quite near zero. For UK all three parameters were better estimated than in ISDW, which is an important result for the aim of the present work.

Another important result is that UK carried out for GRID3 (fewer sampling points) compared with ISDW for GRID1 (higher number of sampling points) showed higher accuracy, fact justified by RMSD and ME values, showing that Kriging, in computational modeling of bathymetric surfaces, is more accurate than ISDW even in unfavorable situations.

In practical terms, one of the reasons for building digital elevation models of water bodies is to subsequently calculate the volumes. Thus, volume calculation of the surveyed reservoir was carried out aiming at verifying the occurrence of significant differences. Results are shown according Table 4.

Table 4:
Calculated volumes for each interpolator and sampling grid

	GRID1	GRID2	GRID3
UK	30,318 m³	30,155 m³	31,118 m³
ISDW	32,689 m³	33,203 m³	33,933 m³

Difference of 2,371 m³ is noticed for GRID1, the difference is even higher for GRID2, around 3,000 m³, yet for GRID3 the difference was approximately 2,800 m³. If we consider UK for GRID1 as the most accurate interpolation, fact justified in Table 3, while applying UK in GRID2 and GRID3 the following discrepancies are found in the volume calculation, respectively: 163 m³ (0.5%) and 800 m³ (2.6%). Whereas while applying ISDW in GRID2 and GRID3, still considering UK is for GRID1 as the most accurate interpolation, the following discrepancies are found in the volume calculation: respectively, 2,885 m³ (9.5%) and 3,615 m³ (11.9%). This is another fact that transmits the control of Geostatistics towards the deterministic method studied here.

Since GRID1 presents more accuracy for both interpolators, the DEMs (Digital Elevation Model) produced from GRID1 for interpolation carried out by UK (a) and ISDW (b) are shown according Figure 7.

Figure 7:
Bathymetric depth DEM based on universal kriging using GRID1 (left map) and Depth bathymetric DEM based on ISDW using GRID1 (right map).

Some differences are noticed by analyzing the DEMs produced by UK and ISDW. The surface produced by universal kriging, for both GRIDS, creates a more uniform floor with smoother outlines. Such result is mainly due to the fact that kriging is an accurate interpolator, which different from ISDW, estimates beyond maximum and minimum limits of sampled point values, without a bias and with minimum variance. It is possible to notice that this interpolator models both regional trends and local anomalies, and is not sensitive to irregularly sampled or grouped data. A disadvantage of this interpolator is the mathematical complexity of its algorithm and necessary time for modeling the semivariogram.

The DEM interpolated by ISDW for both GRIDS presented higher variability when compared with UK. It could also be noticed that ISDW suffers great influence from anomalous local values, besides being sensitive to grouped data, and in their presence, estimates biased values.

4.Conclusions

The construction of bathymetric surfaces is an important component in several studies.

This study allowed verifying that the kriging interpolator presented better results in comparison to the inverse distance interpolator for this set of data; in scattered and abundant sample GRIDS. It was also verified, as it is standard in Surveying Engineering, that the volume calculation of the damming in study was more accurate when UK was applied in a scattered sample GRID, comparatively to the ISDW applied in abundant data. Another reason for the use of kriging is the possibility of generating DEM of uncertainties of the interpolation.

In view of the present results, the application of Geostatistics is recommended in the modeling of bathymetric surfaces, with either scattered or abundant data.

Considering kriging as superior in the construction of bathymetric surfaces, further studies are recommended to define the best sample GRID, in terms of cost and benefit using UK.

References

• Azpurua, M. dos Ramos, K.2010. A comparison of spatial interpolation methods for estimation of average electromagnetic field magnitude. Progress in Electromagnetics Research M, 14, pp.135-145.
• Babak, O. and Deutsch, C. V. 2009. Statistical approach to inverse distance interpolation. Stochastic Environmental Research and Risk Assessment, 23(5), pp.543-553.
• Banerjee, S., Carlin, B. P. and Gelfand, A. E. 2015. Hierarchical Modeling and Analysis for Spatial Data New York: CRC Press.
• Bello-Pineda, J., Stefanoni-Hernández, J. L. 2007. Comparing the performance of two spatial interpolation methods for creating a digital bathymetric model of the Yucatan submerged platform. Pan-American JournalofAquatic Sciences 2, pp.247-254.
• Camargo, E. C. G. 1998. Geoestatística: fundamentos e aplicações In: Camara, G., Medeiros, J.S. Geoprocessamento em projetos ambientais. São José Dos Campos: Inpe.
• Carvalho, J. R. P. and Assad, E. D. 2005. Análise espacial da precipitação pluviométrica no estado de são paulo: comparação de métodos de interpolação. Engenharia Agrícola, 25(2), pp.377-384.
• Chilès, J. P. and Delfiner, P. 2012. Geostatistics: modeling spatial uncertainty 2^nd edition. New York: J. Wiley.
• Curtarelli, M. et al. 2015. Assessment of Spatial Interpolation Methods to Map the Bathymetry of an Amazonian Hydroelectric Reservoir to Aid in Decision Making for Water Management. ISPRS International Journal of Geo-Information 4, pp.220-235.
• Felgueiras, C. A. 1998. Modelagem númerica de terreno In: Camara, G., Medeiros, J.S. Geoprocessamento em projetos ambientais . São José Dos Campos: Inpe.
• Ferreira, I. O., Domingos, R. D. and Santos, G. R. 2015. Coleta, Processamento e Análise de dados batimétricos Novas edições acadêmicas.
• Gotway, C. A. and Hartford, A. H. 1996. Geostatistical methods for incorporating auxiliary information in the prediction of spatial variables. Journal of Agricultural, Biological, and Environmental Statistics, pp.17-39.
• Guimarães, E. C. 2004. Geoestatística básica e aplicada Uberlândia: Núcleo de estudos estatísticos e Biométricos.
• IHO. 2005. Manual on hydrography Mônaco: International Hydrographic.
• Kanegae J. H. et al. 2006. Avaliação de interpoladores estatísticos e determinísticos como instrumento de estratificação de povoamentos clonais de Eucaliptus sp. Revista Cerne, pp.123-136.
• Kravchenko, A. N. 2003. Influence of spatial structure on accuracy of interpolation Methods. Soil Science Society of America Journal, 67(5), pp.1564-1571.
• Landim, P. M. B. 2000. Introdução aos métodos de estimação espacial para confecção de mapas. Revista Laboratório Geomat, pp.02-20.
• Landim, P. M. B., Sturaro, J. R. and Monteiro, R. C. 2000. Krigagem ordinária para situações com tendência regionalizada. Revista Laboratório Geomat, pp.06-12.
• Lark, R. M. 2009. Kriging a soil variable with a simple nonstationary variance model. Journal of Agricultural, Biological and Environmental Statistics, 14(3), pp.301-321.
• Laslett, G.M. et al. 1987. Comparison of several spatial prediction methods for soil pH. Journal Soil Science, 38, pp.325-341.
• Meng, Q., Liu, Z. and Borders, B. E. 2013. Assessment of regression kriging for Spatial interpolation-Comparisons of seven GIS interpolation methods. Cartography and geographic information science, 40, pp.28-39.
• Merwade, V. M., Maidment, D. R. and Goff, J. A. 2006. Anisotropic considerations while interpolating river channel bathymetry. Journal of Hydrology, 331, pp.731-741.
• Morillo Barragán, J. et al. Análisis de calidad de un modelo digital de elevaciones generado con distintas técnicas de interpolación. XIV Congreso Internacional de Ingeniería Gráfica Santander, España - 5-7 junio de 2002.
• Sánchez, J. A. C. 2010 Cartografía Submarina Madrid: Escuela Técnica Superior de Ingenieros Industriales.
• Santos, G. R., Oliveira, M. S,. Louzada, J.M. and Santos, A.M.R.T. 2011. Krigagem simples versus krigagem universal: qual o preditor mais preciso? Energia na Agricultura, 26, pp.49-55.
• Silva, S. A., Lima, J. S. S., Souza, G. S. and Oliveira, R. B. 2008. Avaliação de interpoladores estatísticos e determinísticos na estimativa de atributos do solo em agricultura de precisão. Revista Idesia, 26, pp.75-81.
• Silveira, T. A., Portugal J. L. Sá, L. A. C. M., and Vital, S. R. O. 2014. Análise estatística espacial aplicada a construção de superfícies batimétricas. Geociências, 33(4), pp.596-615.
• Siqueira, D.S., Marques Júnior, J. and Pereira, G. T. 2010. The use of landforms to predict the variabilitu of soil and orange attributes. Geoderma, 155, pp.55-66.
• Souza, E. C. B. 2003. Determinação das variações volumétricas no ISTMO da ilha do mel utilizando PDGPS. Boletim de Ciência Geodésica, 9, pp.53-74.
• Souza, G. S., Lima, J. S. S., Xavier, A. C., Rocha, W. S. D. 2010. Krigagem ordinária e inverso do quadrado da distância aplicados na espacialização de atributos químicos de um argissolo. Revista Scientia Agraria, 11, pp.73-81.
• Tabios, G. Q. and Salas, J. D. 1985. A comparative analysis of techniques for spatial interpolation of precipitation. Water Resources Bulletin, 21, pp.365-380.
• VIEIRA, S. R. 2000. Geoestatística em estudos de variabilidade espacial do solo. In. Novaes, R.F., Alvarez, V.V.H., Schaefer, C.E.G.R. Tópicos em ciências do solo. Revista Brasileira de Ciência do Solo, 1, pp.2-54.
• Warrick, A. W. Z., Hang, R., Harris, M. R. and Myers, D.E. 1988. Direct comparations between kriging and other interpolation validation of flow and transport models for the unsaturated zone. pp.254-326.
• Warrick, A. W., Nielsen, D. R. 1980. Spatial variability of soil physical properties in the field. In: HILLEL, D. Applications of soil physics New York: Academic Press.
• Watson, D. F., Philip, G. M. 1985. A refinement of inverse distance weighted interpolation. Revista Geoprocessamento, 2, pp.315-327.
• Wollenhaupt, N. C., Wolkowski, R. P. and Clayton, M. K. 1994. Mapping soil test phosphorus and potassium for variable-rate fertilizer application. Journal of Production Agriculture, 7, pp.441-448.

Publication Dates

• Publication in this collection
  Sept 2017

History

• Received
  03 Oct 2016
• Accepted
  18 Mar 2017