Spatial Analysis of Pottery Presence at the Former Pobedim Hillfort (an Archeological Site in Slovakia)

Abstract

The aim of this study was a spatial analysis of the pottery occurrence (potsherds) in the acropolis part of the Pobedim hillfort (Slovakia) using two spatial statistical methods (spatial autocorrelation and kriging interpolation) with the help of GIS and their subsequent comparison. To understand the landscape of the study area, seven land use maps were created for different historical years (1783–1785, 1845, 1882, 1956, 1971, 2010 and 2017) confirming that the study area was predominantly utilized as arable land, which was related to advantageous floodplain location between the rivers of Horný Dudváh and Dubová. Using the Moran coefficient of spatial autocorrelation, it was found that there were seven high–high clusters and three high–low clusters representing the occurrence of potsherds. Using the kriging interpolation, three clusters of high concentration were found. Subsequent comparison of both methods revealed three identical areas with high frequency of pottery occurrence indicating places where significant settlement objects were located, such as the dwelling of a wealthy man, pottery workshop and the like. The difference between the areas with high number of potsherds between the two methods is approximately 12%, which indicates an acceptable match between the two methods and their applicability for spatial (geographic)–archaeological research.

1. Introduction

The fortified settlements were built only in a certain period of Slavic culture development. Therefore, in Central Europe the fortified settlements entitled this period as a “hillfort period”. The surroundings of the fortified settlement were created usually by dwellings of peasants who supplied the hillfort elite primarily with the food. The fortified settlements were built either on hills or in a marshy terrain (i.e., lowland forts). Moreover, the fortified settlements had different functions, such as guarding or strategic, or they also played the role of the tribe center, where the tribal chief and his family as well as his armed companions resided. As part of the Christianization, churches, often the stone ones, were later built on the acropolis of the hillforts serving the tribal chief and his relatives. Various craft activities were also found on the periphery of fortified settlements. In connection with these facts, the fortified settlements were quite large, often occupying several hectares.

The studied phenomena have one thing in common: they exist in space, i.e., they are more or less influenced by their immediate but also distant surroundings. They do not exist in isolation from their surroundings. Therefore, when studying phenomena and the dependence between them, it is necessary to look for the factors that have an impact on these phenomena. One such factor is the space and the links associated with it, which is an essential principle of geography. The basic attributes of spatial data include distance, direction and relative location [1]. Spatial data can be divided into three basic types: points, lines and areas.

In order to measure the relationship between phenomena or events separated by certain spatial or temporal sections, spatial statistics uses the spatial autocorrelation. The main idea of spatial autocorrelation is the assertion that if the values of the investigated feature for each pair of areas of a given space are uncorrelated, then there is no spatial autocorrelation of the investigated feature in the system of areas. This argument is based on the so-called Tobler’s first law of geography [2] according to which everything is related to everything else, but closer things are more related than distant things.

The GIS-based interpolation method is an important part of spatial statistics. There are different interpolation methods than can be used to analyze the data. Particularly, the kriging interpolation method can be used to improve the efficiency of spatial data analysis due to its high prediction power. Based on the statistical characteristics of known samples, this method can quantify the spatial autocorrelation between measurement points, improve the data accuracy or highlight the overall distribution trend [3].

In addition, the kriging method uses a semivariogram, i.e., possible combinations between pairs of points, in order to express spatial continuity, which quantifies the correlation weight as a function of distance. In this sense, data, which are closest to a known point, have higher weight during the interpolation [4]. The geostatistical methods, especially semivariogram and kriging method, are described in detail in a number of works, for example, Goovaerts [5], Webster and Oliver [6] and the like.

Many archaeological works dealt with the research of archaeological finds of pottery, for example, Dewar and Marsh [7], Harush et al. [8], Van Valkenburgh et al. [9], Pecci et al. [10] and others, but their findings were based only on archaeological methods. In terms of spatial analysis, archaeological sites are described, for example, in the work of Papworth et al. [11], who dealt with spatial statistics and modeling of archaeological finds using laser scanning methods. In addition, several other studies highlighted recent advances in archaeological predictive modeling using geospatial technologies [12,13,14,15]. Casas and Tema [16] dealt with spatial and temporal analysis of archaeological sites using the European model of the geomagnetic field SCHA.DIF.3K to calculate the uncertainty of dating over the past 3000 years. Furthermore, the work of Niknami et al. [17] described a model of the archaeological landscape in relation to distribution models of archaeological sites using the spatial processes. Other similar works are, for example, Carrer [18], Gibbs [19], Lasaponara and Masini [20], Spurná [21] and so on. Moreover, spatial relationships and distribution of objects are often analyzed by kriging method in other related geoscience disciplines, for example, when analyzing natural hazards [22], landscape changes [23] or retention changes of water reservoirs [24,25].

The aim of this paper is to analyze archaeological findings, particularly quantity of pottery—number of potsherds, from the former Pobedim hillfort (acropolis) in terms of their spatial distribution using spatial statistics methods, namely Moran coefficient of spatial autocorrelation and kriging interpolation, which results were also compared.

2. Research Area

The archaeological site of the Pobedim hillfort is situated in the southwestern part of the Pobedim municipality near the neighboring municipalities of Očkov and Podolie. Pobedim municipality belongs to the Nové Mesto nad Váhom District (NUTS IV) and Trenčín Region (NUTS III). The fortified settlement is located between the rivers of Horný Dudváh and Dubová and belongs to the geomorphological sub-unit of Dolnovážska niva (floodplain) with an elevation around 163 m a. s. l. Along with Bojná, it belongs to two early-medieval sites in Slovakia, which are known for the great occurrence of archaeological finds. The Pobedim hillfort belongs to the lowland forts.

The Pobedim hillfort is already known from the literature in the second half of the 19th century. The first archaeological finds were published in the 1920s and the first sketch of the hillfort is also from this period although the ground plan of the settlement can be seen on the map from the second military mapping (1810–1869) or on the cadastral map from the year 1900. The first archaeological find was found at this place in 1935 and it was represented by 49 iron axe slabs deposited on the quern stone. In 1956, the hillfort was surveyed for the first time in more detail. During the years 1959–1975, the hillfort was not excavated at all [26].

According to Bialeková [27], the Pobedim hillfort was built at the end of the 8th century and destroyed after a fire at the turn of the first and second thirds of the 9th century, which is probably related to the joining of the Principality of Nitra and Principality of Moravia.

The Pobedim hillfort can be divided into two parts (acropolis and nearby settlement), which probably were not built at the same time. The acropolis covers an area of approximately 4.1 ha and the nearby settlement has an area of approximately 3.9 ha. The geophysical survey revealed also some other enclosed formations connected to the fort system, but their dating is unknown and it is assumed that they were represented by cattle pens and the like [26].

Most of the buildings were found to be close to the rampart while southern part was less inhabited. In addition to iron objects, the most common finding is pottery, quern stones for grinding grain, grinding stones and the like. It is interesting that a large number of quern stones were brought to this place from the quarries in the so-called Slovak Gate in the Pohronie region. Moreover, numerous conical clay cones document the textile production. One of the largest collections of pottery (about 12,000 potsherds) in Slovakia from the 9th century comes from this place (kitchenware, massive grain roasters and the like) [26].

This paper focused, especially, on the acropolis part of the hillfort. Altogether, six sites (I to VI) were analyzed in the acropolis part. Each site was divided into sectors (mostly squares with a side of 5 m) where archaeological excavations were carried out (Figure 1). The number of analyzed potsherds in each sector, based on the archaeological research, is shown in Figure 2.

3. Methods

In the first step, the maps of land use from different years were created to map the spatio-temporal changes of the study area. Spatio-temporal analysis of landscape structure is a common method of assessing land use changes [28,29]. Using ArcGIS 10.2 (Esri, Redlands, CA, USA) software and manual digitization, the following historical military and topographic maps as well as orthophotos were interpreted:

• (a). Map from the first military survey (years 1783–1785) at a scale 1:28,800,

• (b). Map from the second military survey (1845) at a scale of 1:28,800,

• (c). Reambulated map (1936) based on the map from the third military survey (1882) at a scale of 1:25,000,

• (d). Topographic map from 1956 at a scale of 1:25,000,

• (e). Topographic map from 1971 at a scale of 1:25,000,

• (f). Orthophotos from 2010 at a scale of 1:2000 and

• (g). Orthophotos from 2017 at a scale of 1:2000.

3.1. Spatial Autocorrelation

After the analysis of land use maps, the archeological finds of potsherds from the acropolis part of the Pobedim hillfort were analyzed in terms of their spatial distribution. In particular, the Moran coefficient of spatial autocorrelation was used to assess the occurrence of potsherds.

Spatial autocorrelation has a significant role in the study of spatial statistics or spatial analysis [30] and it is defined as the presence of the spatial structure of mapped variables with respect to their geographical proximity [31]. Spatial autocorrelation is a specific type of correlation where the relation of one variable in space and time is evaluated within one observation [32]. There are currently several measures of spatial autocorrelation. They can be divided into the following three groups: global indicators of spatial association, local indicators and variogram-based geo-statistical tools [33].

The model of spatial association (autocorrelation) is given by the following Equation (1) [1]:

(1)$X_i=\rho {\displaystyle \sum }_{j=1}^n{w}_{ij}{X}_j+{\epsilon }_i,i=1, 2,\dots , n,$

where ε_i (i = 1, 2, …, n) are independent equally distributed random variables with common scattering σ²; w_ij (i, j = 1, 2, …, n) are weights specifying the relationship of area i and j; ρ is the measure of the overall level of spatial autocorrelation X_i and X_j for w_ij > 0.

If the value ρ = 0, then there is no spatial autocorrelation. If ρ > 0, we speak about a positive spatial autocorrelation while if ρ < 0, it is a negative spatial autocorrelation.

The positive spatial autocorrelation means that geographically close values of the studied variable are grouped side by side in space and thus tend to group with similar values of the studied variable on the map. To be more concrete, high values are located near high values, medium values near medium values and low values near low values [34]. The negative spatial autocorrelation means that different values of the studied variable are grouped side by side. The absence of spatial autocorrelation means that the values of the studied phenomenon occur randomly in space. In this study, the null hypothesis H₀ was tested: there is no spatial autocorrelation (distribution of the number of archaeological finds of pottery in the acropolis is random) against the alternative hypothesis H₁: there is a spatial autocorrelation (distribution of the number of archaeological finds of pottery in the acropolis is not random).

When testing the existence of spatial autocorrelation, the Moran indicator of spatial autocorrelation is mostly applied using the Equation (2):

(2)$I=\frac{n}{S_0}·\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{w}_{ij}\left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}{{\sum }_{i=1}^n{w}_{ij}{(x_i-\overline{x})}^2},S_0={\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=1}^n{w}_{ij},$

where n is the number of spatial units (number of sectors); w_ij is spatial weight while w_ij = 1, if sectors i and j are neighboring, w_ij = 0 otherwise (i, j = 1, 2, …, n); x_i (i = 1, 2, …, n) is the value of the studied phenomenon (quantity or frequency of potsherds) in the sector I and $\overline{x}$ is the arithmetic mean of the studied phenomenon given by the Equation (3):

(3)$\overline{x}=\frac{{\sum }_{i=1}^n{x}_i}{n}.$

The Moran coefficient of spatial autocorrelation I takes values from −1 to +1 [35]. The closer the value of the Moran index I is to 0, the more randomness is indicated, i.e., statistical insignificance of the variable in space. The closer the value of I is to 1, the more positive spatial autocorrelation is indicated and if I is closer to −1, the negative spatial autocorrelation is indicated.

Statistical significance of the Moran index of spatial autocorrelation needs to be verified [36]. One of the significance tests, which verifies the statistical significance of the Moran coefficient of autocorrelation, i.e., the positive or negative spatial autocorrelation, is based on a comparison of the z-score value and the corresponding probability value of p. The z-score is given by Equations (4)–(6):

(4)$Z_I=\frac{I-E\left(I\right)}{\sqrt{V\left(I\right)}},$

where:

(5)$E\left(I\right)=\frac{-1}{n-1},$

(6) $V\left(I\right)=E\left(I^2\right)-E{(I)}^2.$

In practice, the greater the z-score value, i.e., the more distant it is from the value of 0, the greater the level of confidence that the studied phenomenon is autocorrelated.

In addition to the z-score value, the corresponding value of p is also computed using ArcGIS 10.2 software, i.e., the probability with which the analyzed phenomenon is a result of the random process. The p value is thus the probability of an error, which is made by rejecting the tested hypothesis. Therefore, the test results can be evaluated based on the calculated probability values of p. If the calculated p value is sufficiently small (p < 0.05 or p < 0.01), the tested hypothesis (Moran coefficient of spatial autocorrelation is not statistically significant) at a significance level of 0.05 or 0.01 can be rejected. Otherwise, the tested hypothesis cannot be rejected [35,36].

Moreover, it should be noted that the calculation of global characteristics of spatial autocorrelation is based on the assumption of homogeneity. If this assumption is not met, global statistics may misrepresent the absence of spatial autocorrelation in the analyzed dataset although, in fact, there is a strong positive autocorrelation in one part of the locality and a strong negative autocorrelation in another part of the locality. In view of the above reasons, it is appropriate to use the local indicators of spatial association, i.e., the so-called LISA (Local Indicators of Spatial Association) that are site-specific. Anselin [37] defined the LISA as a statistics that fulfills the following two conditions:

• (1). The LISA for each area (observation) indicates the extent of significant spatial clustering of similar values around this area.

• (2). The sum of LISA for all observations is proportional to the global indicator of spatial association.

According to Anselin [38], five different options may occur in the studied area when using the LISA:

• (1). Places with high values and similar neighbors: (high–high or H–H), known as hot spots, showing a scenario of positive spatial autocorrelation,

• (2). Places with low values and similar neighbors: (low–low or L–L), called cold spots, showing also a scenario of positive spatial autocorrelation,

• (3). Places with high values and neighbors with low values: (high–low or H–L), called potential spatial outliers, showing a negative spatial autocorrelation,

• (4). Places with low values and neighbors with high values: (low–high or L–H), again called spatial outliers, showing a negative spatial autocorrelation,

• (5). Places with no significant local spatial autocorrelation (not significant).

The LISA analysis is closely related to the Moran plot through which the basic results of spatial autocorrelation analysis can be shown. Anselin [38] was the first who introduced the so-called Moran scatter plot with points and coordinates $(x_i,\sum _j{w}_{ij}{x}_j)$, which allows us to reveal significant local structures. In addition, the slope of these points is the Moran coefficient of spatial autocorrelation. Since the variables are taken as deviations from their averages, the Moran diagram is centered at the position 0.0. Four quadrants in the Moran diagram represent four different types of relationships between the original values of the variable located on the horizontal axis and the average values of the neighboring units located on the vertical axis. The quadrants in the Moran diagram indicate how much each type of spatial dependence impacts the generation of the stated value of global Moran statistics. Units located in the top right and bottom left quadrants represent positive spatial associations. In case of the upper right quadrant, we speak about the relation between high values (hot spots) while in the lower left quadrant, it is the relation between low values (cold spots). Negative associations correspond to the upper left and lower right quadrants, i.e., low values surrounded by high values (upper left) and high values surrounded by low values (bottom right)—spatial outliers [38]. The Moran scatter plot is also commonly used to test the statistical significance of local coefficients of spatial autocorrelation. In this study, the Moran scatter plot was created in GeoDa software (The University of Chicago, Chicago, IL, USA).

3.2. Kriging Interpolation

Another method used for analyzing the spatial distribution of pottery presence in the acropolis of the Pobedim hillfort was the kriging interpolation performed in ArcGIS 10.2 software.

The kriging interpolation is an advanced geo-statistical method, which is based on statistical models, which contain autocorrelation, i.e., statistical relationships between the measured points. On one hand, geo-statistical methods have the ability to produce a prediction surface while on the other hand, they provide some degree of certainty or accuracy of performed predictions [39].

According to Burrough [39] and Oliver [40], the kriging method is based on the precondition that the distance or direction among the studied points reflects a spatial correlation, which is then used to explain the change in surface. The kriging interpolation performed in ArcGIS software adjusts the mathematical function to a specified number of points or to all points within a specified radius in order to determine the resulting value for each location. The GIS-based kriging interpolation process consists of the following tasks: (a) includes exploratory statistical data analysis; (b) variogram modeling and (c) creating the surface and (where appropriate) investigating the surface variance.

The kriging method weighs surrounding measured values to derive a prediction for an unmeasured location. The general formula (Equation (7)) for the kriging method is as follows [39,40]:

(7)$\widehat{Z}(s_0)={{\displaystyle \sum }}_{i=1}^N{\lambda }_{i}Z(s_i),$

where Z (s_i) is the measured value of the i-th location, λ_i is the unknown weight for the measured value of the i-th location, s₀ is the prediction location and N is the number of measured values.

Based on Equation (7), weights are related to the distance between the measured points as well as the prediction location and the overall spatial arrangement of the measured points. Therefore, spatial autocorrelation have to be quantified in order to use the spatial arrangement of the weights. As for the ordinary kriging interpolation, the weight λ_i depends on the model, which is applied to the measured points, distance to the prediction location as well as spatial relationships among the measured values in radius of the prediction location [39,41].

When applying the spatial modeling of the structure of the measured points using the kriging method, a graph of an empirical semivariogram is usually constructed. The following Equation (8) is used to calculate the semivariogram [42]:

Semivariogram (distance_h) = 0.5 × average ((value_i − value_j)²).(8)

According to McBratney and Webster [42], the use of semivariogram is an inevitable task, which is performed between spatial description and spatial prediction. The kriging interpolation has the advantage to predict attribute values at undetected locations while the empirical semivariogram provides information on spatial autocorrelation of datasets. The kriging interpolation available in ArcGIS software enables to choose from the following functions for modeling the empirical semivariogram: spherical, circular, exponential, Gaussian and linear. The prediction of unknown values is influenced by the selection of one of these functions.

The kriging interpolation differs between the ordinary and universal kriging. The ordinary kriging, which was also used in this study, is the most common kriging technique. It assumes that the constant mean is not known. According to Manjarrez-Domínguez et al. [4], the formula that expresses the ordinary kriging interpolation method is as follows (Equation (9)):

(9)${\displaystyle \sum }_{i=1}^n{\lambda }_iƔ\left[d\left(s_i,s_j\right)\right]+m=Ɣ\left[d\left(s_i,s_j\right)\right]$

where n is the number of observations, m is the Lagrange multiplier that is used to minimize constraints, λ_i is the weight assigned to each observation (the sum of all weights equals 1), i, j are observations, 0 is the point of estimation, s is the point of estimation and d(s_i, s_j) is the distance between s_i and s₀ from the semivariogram.

4. Results

After creating the land use maps for the studied years, it was found that the Pobedim archaeological site was utilized mainly as agricultural land (fields; Figure 3). The reason for agricultural use of this site is mainly because of fertile soils, namely mollic fluvisols and mollic gleysols, which are heavy (clay–loam) and very heavy (clayey and clay) soils. Agricultural land is interrupted only by a tree-line. Based on the maps of first and second military surveys, a small area of grassland was also interpreted. Since the maps from these periods are not that accurate, especially the map from the first military survey, which was constructed without any astronomical-geodetic measurements, the grasslands did not have to be included in the study area at all. In 1959, the Pobedim archaeological site was destroyed by an agricultural cooperative as part of land consolidation, which can be seen in Figure 3e–g.

Furthermore, the Moran coefficient of spatial autocorrelation was implemented in ArcGIS 10.2 software in order to assess the spatial distribution of potsherds found in the acropolis. The result of the spatial autocorrelation is presented in Table 1 and Figure 4.

The null hypothesis was tested, i.e., hypothesis (H0): there is no spatial autocorrelation of the number of archaeological finds of pottery in the acropolis or their spatial distribution is random.

On the contrary, the alternative hypothesis was hypothesis (H1): there is a spatial autocorrelation of the number of potsherds in the acropolis or their spatial distribution is not random.

The testing of validity of the null hypothesis was also based on the local Moran coefficient of spatial autocorrelation.

Based on the results of the Moran coefficient of spatial correlation, which is shown in Table 1, it can be concluded that the calculated probability value of p is p = 0.0000. Since the p-value was in all cases less than 0.01, the null hypothesis was rejected at the significance level equal to 0.01 and the alternative hypothesis H1 was accepted, i.e., the spatial distribution of potsherds in sectors in the acropolis is not random.

In addition to the z-score, the statistical significance of the Moran coefficient was tested also using the Moran scatter plot (Figure 5), which was created in the GeoDa software.

As already mentioned in Section 3.1., the quadrants of the Moran diagram indicate the share of individual types of spatial dependence on the generation of the calculated value of Moran statistics. It is apparent from Figure 5 that there was a significant correlation between the number of potsherds in the sector and the number of potsherds in the neighboring sectors. At the same time, it could be claimed that, to a large extent, there were spatially grouped high values, i.e., large numbers of potsherds (high–high quadrant, H–H), as well as spatially grouped low values, i.e., low numbers of potsherds (low–low quadrant, L–L). The membership of the high–low quadrant (H–L) and the low–high quadrant (L–H) indicates spatial outliers, i.e., sectors in which significantly higher or lower number of potsherds was recorded compared to the neighboring sectors.

Based on the Moran scatter plot as well as the calculated corresponding z-score values, it can be concluded that the Moran coefficient of spatial autocorrelation was statistically significant for the Pobedim acropolis. The statistical conclusiveness of the Moran coefficients of spatial autocorrelation means that there are sectors next to each other that are closer to each other in terms of the number of potsherds. Therefore, the number of potsherds found in a given sector is dependent on the number of potsherds in the surrounding sectors. This means that the frequency of archaeological finds is influenced by the location of the sectors.

Since the Moran coefficient of correlation was statistically significant for the number of potsherds, the next step was to investigate which sectors in the Pobedim acropolis cause spatial autocorrelation for the abovementioned frequencies. The statistically proven local Moran coefficients for the number of potsherds in the Pobedim acropolis are shown for individual finds (sectors) in Figure 6.

Based on Figure 6, it could be seen that the local Moran coefficient identified seven clusters of sectors with high number of potsherds, i.e., there was a positive spatial autocorrelation of a high number of potsherds in the acropolis. These twenty-five high–high sectors were located in sites I, II, III and IV. The site II contained three clusters: the first largest cluster of nine sectors was located in the western part, the second cluster with two sectors in the central part and the third cluster with one sector in the northeastern part, also forming one complex cluster with the cluster, which was identified in the western part of the site I (two sectors). Furthermore, in site I there were two clusters of two sectors in the central part. Site III had one cluster with three sectors located in the northwestern part. Site IV had one cluster with four sectors located in the northern part.

Moreover, the results of the LISA analysis point to three high–low sectors, which was the negative spatial autocorrelation. These sectors had a high number of potsherds that were adjacent to sectors with low number of potsherds. All high–low clusters were located on site I.

It could be seen that the local Moran coefficient of spatial autocorrelation of the number of potsherds in other sectors of the acropolis was not statistically significant, i.e., these sectors did not show significant spatial autocorrelation. The potsherds found at these sectors were randomly arranged having different numbers.

After performing the spatial analysis using the kriging interpolation method, three areas with high concentration of potsherds were found in the acropolis of Pobedim hillfort (Figure 7). The first cluster was located in the eastern part of site I, the second cluster was found in the northwestern part of the site III and the last cluster was located in the northwestern part of the site IV. The spatial analysis performed by the kriging interpolation method did not record a high concentration of potsherds on site I, as compared to the spatial autocorrelation method.

Based on the comparison of the map of spatial autocorrelation and the map of spatial analysis by the kriging interpolation, almost identical areas with high number of potsherds in the acropolis of Pobedim hillfort were identified (Figure 8). As for the map of spatial autocorrelation, the clusters with a high number of potsherds (high–high) cover an area of 625 m² while these clusters in the map created by the kriging interpolation had an area of almost 551 m². The difference between the two methods was less than 12%. These findings indicate that in these localities (site II—western part, site III—northwestern part and site IV—northwestern part), there were important historical settlements, such as the dwelling of a wealthy man, pottery workshop and the like. However, the most important finding, in view of the archaeological research in this locality, was the largest cluster of sectors (spatial autocorrelation: 225 m²) as well as high and very high concentration—last four intervals (kriging: 194 m²; Figure 6, Figure 7 and Figure 8) located in the western part of site II. The reason was that during the fire of the hillfort in these sectors, the remains of fortifications slid to the dwelling and thus remained in-situ. The significance of this archeological finding was confirmed by both research methods applied in GIS: spatial autocorrelation and kriging interpolation.

5. Discussion

As stated by Niknami et al. [17], already people in prehistoric times built their dwellings in slightly elevated areas or near watercourses. This fact is also valid for the selected study area, where people started to build their dwellings at the turn of the 8th and 9th centuries choosing the slightly elevated and fertile area between two rivers. These authors also state that the spatial grouping of archaeological sites is determined based on two reasons. It is namely the nature of the landscape and cultural and socio-economic factors. The nature of the landscape, i.e., elevation, slope, landscape structure, availability and persistence of water resources and food as well as the quality and quantity of transport networks, also influenced the construction of dwellings in the study area. As mentioned in Section 2, a large number of stone quarries were found at the selected archeological site, which was brought here from quarries in the Slovak Gate at the Pohronie Region. The Slovak Gate is approximately 90 km far from the Pobedim archaeological site. This finding suggests that the inhabitants of the Pobedim hillfort were mobile and that there had to be a good transport connection between these two places. On the other hand, cultural and socio-economic factors help to explain the overall nature of archaeological sites and the strong link between the locations of sites [43].

To correctly identify the nature of the landscape in the study area, the historical maps from the first, second and third military surveys, topographic maps from the years 1956, 1971 and also orthophotos from 2010 and 2017 were analyzed. This analysis confirmed that the study area was economically utilized almost during the whole studied period. The advantages of analyzing the historical, topographic maps and orthophotos for the research of archaeological findings were also pointed out in the works of Zupancich et al. [44], Gennaro et al. [45], Tache et al. [46], Drap et al. [47] and so on.

The boundaries of the Pobedim hillfort were already marked on the maps of the second military survey. However, these boundaries were not marked in the map of the first military survey. This can be explained by the fact that maps from this survey did not undergo any astronomical-geodetic measurements. Moreover, the absence of the mathematical bases of these maps resulted in incorrect geometry of the objects and the inaccurate position of objects reduced the informative value and usability of this map [48,49].

Placing the boundary of the study area is another critical research problem. The same points appear regular in the context of the study area, i.e., evenly distributed. On the other hand, in the context of differently delimited area, points were clustered. It follows that no measure of variability should be interpreted independently of the study area [1].

In this study, the spatial statistical methods were used to analyze the spatial distribution of archaeological finds at the former Pobedim hillfort. The spatial autocorrelation method revealed seven high–high clusters, which were described in Section 4. Subsequently, the spatial analysis with the kriging interpolation method was applied. The results of these two methods were very similar suggesting that they were correctly applied for assessing the spatial distribution of potsherds in the study area. The methods of spatial autocorrelation and kriging interpolation were also used in many previous studies, such as [50,51,52,53]. However, it should be noted that these works did not deal with archaeologically based research and they were focused on geographic, environmental or ecological research. Based on the review of literature, similar research to the one presented in this study using selected methods of spatial analysis has been carried out only in few works so far. For instance, Nakoinz and Knitter [54] dealt with the integration of geography principles and archaeology in terms of modeling human behavior in landscapes. They pose different questions related, for example, to the choice or application of spatial-statistical (modeling) methods, which is not arbitrary and should be based on the research objectives, data and theory. According to these authors, the research question, which is studied by two methods, does not allow a specific interpretation, especially, when the two methods produce different results. However, in the presented research the applied methods of spatial autocorrelation and kriging interpolation resulted in very similar results (less than 12% difference), which enable specific interpretations of the results.

The similarity between the results of the spatial autocorrelation and kriging interpolation lies predominantly in the fact that kriging is inherently based on spatial autocorrelation characteristics [55]. Furthermore, the quantification of the spatial autocorrelation’s strength and reach in the observed data was calculated with the semivariogram. According to Olea [56] and Wackernagel [57], kriging interpolation through variography provides an optimal interpolation estimate from observed values and their spatial relationships. Nakoinz and Knitter [54] further point out that due to the semivariogram, it is not necessary to apply a test of spatial autocorrelation.

In addition, the similarity between the two methods is seen in the fact that both methods recognize that data, which are close together are more likely to be similar than data that are farther apart, which is the Tobler’s first law of geography [2]. On the other hand, the main difference between the two presented methods is that spatial autocorrelation was performed for individual sectors (mostly squares where the archeological excavations were carried out) while the result of the kriging interpolation is the continuous surface of values, which is caused by the fact that kriging is a contouring method.

6. Conclusions

The study area was mainly agriculturally utilized throughout the studied period (from 1783 to 2017). Due to the agricultural purposes, the archaeological site was destroyed in 1959.

As mentioned in Section 4, spatial autocorrelation defined seven clusters. These clusters indicate places where significant settlement objects were located, such as the dwelling of a wealthy man, pottery workshop and the like. Especially, the cluster in the western part of site II is very important because archaeological research has revealed that during the fire of the Pobedim hillfort, the remains of the fortification slid to the dwelling and thus the finds were preserved in-situ. As a result, eight reconstructable pots and several iron objects (iron spurs, agricultural and craft tools) were found in this object. It was interesting that some of the iron objects were well-preserved and those that were damaged could have represented a raw material for the production of these iron objects. It is noteworthy that six pots and iron objects were found in such a position as if they had originally been placed on a hanging shelf and had fallen to the ground due to the fire of the object. Other high–high clusters were also found in the site III. This cluster with high number of potsherds could have been related to the localization of the production facility or dwelling in the acropolis. Regarding the statistically significant high number of potsherds in site IV, their occurrence is random from the archaeological point of view.

A total of twenty one settlement objects were uncovered in the acropolis part of the Pobedim hillfort. Some of the objects served as dwellings equipped with a stone furnace while some of them were production objects and the like. The total number of potsherds found at Pobedim hillfort was 12,000 pieces. After the cessation of the hillfort, its area was used for skeletal burial. It is assumed that during the fire that destroyed the hillfort, the walls slid to the dwellings and due to that many potsherds were preserved in-situ. Since the archaeological research could not continue because of the disagreement of the local agricultural cooperative, it remains questionable what was placed in the central part of the acropolis.

Author Contributions

Conceptualization, J.V., A.T., M.V. and I.V.; methodology, J.V., A.T. and M.V.; software, J.V. and M.V.; formal analysis, J.V. and A.T.; resources, J.V. and I.V.; data curation, J.V. and M.V.; writing—original draft preparation, J.V., A.T., M.V. and I.V.; writing—review and editing, M.V. and J.V.; visualization, J.V. and M.V. All authors discussed the results and edited the manuscript.

Funding

This work was supported by the Slovak Research and Development Agency under the Contract no. APVV-18-0185 and by the VEGA agency (Scientific Grant Agency of the Ministry of Education, Science, Research and Sport of the Slovak Republic and the Slovak Academy of Sciences) under the grant no. 1/0934/17.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Table

Generate image description

Figure 1. Research area of the Pobedim hillfort.

Generate image description

Figure 2. Number of potsherds in the acropolis part of the Pobedim hillfort.

Generate image description

Figure 3. Land use: (a) first military survey, (b) second military survey, (c) third military survey (reambulated map), (d) topographic map 1956, (e) topographic map 1971, (f) orthophoto 2010 and (g) orthophoto 2017.

Generate image description

Figure 4. Result of the Moran coefficient of spatial autocorrelation, z-score and p-value using ArcGIS 10.2 software.

Generate image description

Figure 5. Moran scatter plot for the number of potsherds in the acropolis (Pobedim hillfort).

Generate image description

Figure 6. Statistically proven local Moran coefficients for the number of potsherds in individual sectors of the acropolis.

Generate image description

Figure 7. Spatial distribution of potsherds in the acropolis of the Pobedim hillfort using the kriging interpolation.

Generate image description

Figure 8. Comparison of the results of spatial autocorrelation and kriging interpolation.

Table 1

Result of testing the spatial autocorrelation of the number of potsherds in the acropolis (Pobedim hillfort).

Variable	Moran Coefficient	p-Value	z-Score Value
Potsherds (number)	0.4410	0.0000 *	9.3545

* Test is statistically significant.