Evaluation of Landsat 8 OLI/TIRS Level-2 and Sentinel 2 Level-1C Fusion Techniques Intended for Image Segmentation of Archaeological Landscapes and Proxies

The use of medium resolution, open access, and freely distributed satellite images, such as those of Landsat, is still understudied in the domain of archaeological research, mainly due to restrictions of spatial resolution. This investigation aims to showcase how the synergistic use of Landsat and Sentinel optical sensors can efficiently support archaeological research through object-based image analysis (OBIA), a relatively new scientific trend, as highlighted in the relevant literature, in the domain of remote sensing archaeology. Initially, the fusion of a 30 m spatial resolution Landsat 8 OLI/TIRS Level-2 and a 10 m spatial resolution Sentinel 2 Level-1C optical images, over the archaeological site of “Nea Paphos” in Cyprus, are evaluated in order to improve the spatial resolution of the Landsat image. At this step, various known fusion models are implemented and evaluated, namely Gram–Schmidt, Brovey, principal component analysis (PCA), and hue-saturation-value (HSV) algorithms. In addition, all four 10 m available spectral bands of the Sentinel 2 sensor, namely the blue, green, red, and near-infrared bands (Bands 2 to 4 and Band 8, respectively) were assessed for each of the different fusion models. On the basis of these findings, the next step of the study, focused on the image segmentation process, through the evaluation of different scale factors. The segmentation process is an important step moving from pixel-based to object-based image analysis. The overall results show that the Gram–Schmidt fusion method based on the near-infrared band of the Sentinel 2 (Band 8) at a range of scale factor segmentation to 70 are the optimum parameters for the detection of standing visible monuments, monitoring excavated areas, and detecting buried archaeological remains, without any significant spectral distortion of the original Landsat image. The new 10 m fused Landsat 8 image provides further spatial details of the archaeological site and depicts, through the segmentation process, important details within the landscape under examination.


Introduction
Object-based image analysis (OBIA) builds on widely used remote sensing techniques such as image segmentation, edge detection, feature extraction, and classification concepts, having been originally developed as a result of the recognition of limitations working with pixel-based image approaches [1,2]. While OBIA has been thoroughly studied in the past concerning multispectral satellite image segmentation [3][4][5][6][7], this concept has been only recently introduced in the domain of archaeological research [8]. As argued by [9], OBIA is a useful and potential tool for archaeological research; however, it is essential to support the development of approaches that integrate OBIA with The paper is organized as follow: The next section presents the archaeological site selected as a case study, the satellite images used, as well as the various algorithms and equations for the fusion modeling and the segmentation analysis; the results are presented in Section 3; while the paper ends with some general conclusions and future work.

Area of Interest
The archaeological site of "Nea Paphos", located at the western part of Cyprus, was selected as the case study. The site is enlisted as a World Heritage Monument, due to its significant archaeological findings and materials. According to written sources, the city of "Nea Paphos" was founded at the end of the 4th century while in the beginning of the next century, when the island of Cyprus became part of the Ptolemaic kingdom. "Nea Paphos" became the center of Ptolemaic administration on the island. Until the end of the 2nd century B.C., "Nea Paphos" was considered to be an important political and economic center of the region [33].
The selection of this site was based, to a large extent, on its ablitity to be detected from both Landsat and Sentinel sensors. The site of "Nea Paphos" includes excavated areas, as well as interesting archaeological proxies, indicators of potential underground archaeological remains which are not yet excavated (see [34] and Figure 1b, white arrows). The site is surrounded in the western and southern part by the Mediterranean Sea, while on the northern and eastern part by the modern city of Paphos. The area is depicted in Figure 1, in a high-resolution aerial orthophoto image.

Datasets
For the purposes of the study, a recent Landsat 8 OLI/TIRS scene (Scene Id: LC81770362019351LGN00) and a Sentinel 2 image (Granule Id: L1C_T36SVD_A023470_20191220T083342) were downloaded from the USGS Earth Explorer platform. Both images were set with minimum cloud coverage (less than 10%) while the time of acquisition between these two images differed by only three days (17 December 2019 for the Landsat

Datasets
For the purposes of the study, a recent Landsat 8 OLI/TIRS scene (Scene Id: LC8177 0362019351LGN00) and a Sentinel 2 image (Granule Id: L1C_T36SVD_A023470_20191220T083342) were downloaded from the USGS Earth Explorer platform. Both images were set with minimum cloud Remote Sens. 2020, 12, 579 4 of 19 coverage (less than 10%) while the time of acquisition between these two images differed by only three days (17 December 2019 for the Landsat image and 20 December 2019 for the Sentinel image). The zenith and sun azimuth angles of both images were comparable, thus, providing analogous solar illumination (sun azimuth~160 • ), and sensor viewing angle (zenith angle~60 • ). The Landsat 8 image was downloaded at a Level 2 processing, which corrected the image both in terms of geometry and radiometry, including corrections for atmospheric effects. The Sentinel 2 image was downloaded at the Level 1C processing providing top-of-atmosphere (TOA) reflectance data.
The USGS Earth Explorer provides corrected images (Level 2) of Landsat 8 for the first seven spectral bands (Band 1 to Band 7) which cover the coastal aerosol, blue, green, red, near-infrared (NIR), and short-wavelength infrared (SWIR-1 and SWIR-2) part of the spectrum, respectively. The spatial resolution of the above-mentioned spectral bands is 30 m. In contrast, for the fusion needs, the Sentinel 2 Level 1C images were restricted to only the available 10 m spatial resolution bands which include Band 2 to Band 4 and Band 8, and correspond to the blue, green, red and near-infrared part of the spectrum.
It should be mentioned that the images at the specific level of processing (i.e., Level 2 for Landsat and Level 1C for Sentinel images) are currently provided as "ready products", meaning that these series of datasets can be directly retrieved from earth data engines such as those of the USGS Earth Explorer platform, without any preprocessing from the end users. Although Sentinel 2 Level 2A (bottom of atmosphere (BOA) reflectance and atmospheric corrected) images are currently not available and not provided through the USGS Earth Explorer platform, such images can be retrieved from the Sentinel Hub engine [35]. Nevertheless, these images are only retrieved at a 20 m spatial resolution [36], which is not efficient in our case study.
In our case study, the specific image (Sentinel 2 Level 1C) was acquired on a clear non-hazy day with minimum cloud cover, thus, minimizing potential impact from the atmosphere. The contribution of the atmosphere was less than 0.2% (reflectance), based on the darkest pixel (DP) image-based atmospheric correction algorithm. In direct comparison with the Sentinel 2 Level 2A (of the same date, downloaded from Sentinel Hub at 20 m pixel resolution), the difference between corrected and non-corrected atmospheric Sentinel images for Bands 2, 3, 4 and 8 was estimated at 0.03% for Bands 2, 3, and 8 and 0.05% for Band 4, which was considered not important.

Methodology
Using the four 10 m spectral bands of the Sentinel 2 image, as the high-resolution input image, fusion pansharpening techniques were implemented to improve the spatial resolution of the Landsat image. At this step, four different pansharpening methods, namely the Gram-Schmidt, the Brovey, the principal component analysis (PCA), and the hue-saturation-value (HSV) algorithms, were implemented.
The Gram-Schmidt method is based on a general algorithm for vector orthogonalization, i.e., the Gram-Schmidt orthogonalization. The Gram-Schmidt fusion simulates the high-resolution band from the lower spatial resolution spectral bands. In general, this is achieved by averaging the multispectral bands. As the next step, a Gram-Schmidt transformation is performed for the simulated high-resolution band and the multispectral bands, where the simulated high-resolution band is employed as the first band [37,38].
Brovey transformation is a simple but widely used red, green, and blue (RGB) color fusion [39]. Brovey transform is one of the most commonly used pansharpening methods due to its high degree of spatial enhancement, speed, and ease of implementation [40]. Since this transform is intended to produce RGB images, only three spectral bands at a time can be merged from the multispectral input scene. The Brovey transform is shown in Equation where H, R, G, and B are for high-resolution, red, green, and blue bands, respectively.
Principal component analysis (PCA) is a statistically-based transformation [41]. The PCA converts intercorrelated multispectral bands into a set of uncorrelated components using orthogonal reprojections of the spectral space. The high-resolution image is fused into the low-resolution multispectral bands by performing a reverse PCA transform [39,41]. The method can be simplified as followed [42]: where F i is the i-th band in the fused (pansharpened) image, MS i ↑ is the upscaled i-th multispectral (MS) band, vi represents the i-th element in the most significant eigenvector, H is the high-resolution image, and PC1 represents the first principal component in the transformation of the multispectral image.
In the hue-saturation-value (HSV) fusion model, hue (H) defines pure color in terms of "green", "red", or "magenta", while saturation (S) defines a range from pure color (100%) to grey (0%) at a constant lightness level. Finally, value (V) refers to the brightness of the color. The HSV method transforms the RGB image to the HSV color space, while it replaces the value band with the high-resolution image. Then, the method resamples the hue and saturation bands to the high-resolution pixel size and transforms the image back to the RGB color space. More details for all the above fusion methods are found at [42,43]. The implementation of all these fusion methods was carried out in the ENVI Harris Geospatial Solutions v.5.1 environment.
The outcomes of the fusion between the Landsat and the Sentinel images were evaluated with various quality image methods. In this study five different methods were used, namely the (a) bias, (b) image entropy, (c) ERGAS (erreur relative globale adimensionnelle de synthèse), (d) RASE (relative average spectral error), and (e) RMSE (root mean squared error). Their equations are given below, i.e., Equations (3) to (7). The bias method provides the deviation degree between the fused and the low-resolution images, while image entropy quantifies the information of the fused image based on Shanon theorem. The ERGAS estimates the ratio between pixels of the fused and the low-resolution image, the RASE method characterizes the average performance of the image fusion for all spectral bands considered, and the RMSE measures the difference between the reference image and the fused image. More details regarding the aforementioned quality image methods are found in [44,45]. The implementation of the quality methods was carried out in the MathWorks MATLAB R2016b environment based on [45] toolbox.
On the basis of the above findings, the best fusion model and spectral 10 m band were selected. The next step included the segmentation process in the ENVI Harris Geospatial Solutions v.5.1 environment. Several segments using the edge method were implemented with various scale factors. The latest included segmentations from a scale factor of 10 to 100, with a step of 10. The segmentation process (edge method) produces a gradient image using the Sobel edge detection method. Then, the watershed algorithm is applied to the gradient image [46]. The watershed segmentation is a Remote Sens. 2020, 12, 579 6 of 19 region-based method that has its origins in mathematical morphology [47]. In watershed segmentation, an image is regarded as a topographic landscape with ridges and valleys [48].
The segmentation performance, as a comparison of the segments produced from the 30 m Landsat 8 image and the fused 10 m Landsat 8 image was then examined. The first image was used as the ground truth image, while the second image was used as the slave image. This procedure was carried out once again in the MathWorks MATLAB R2016b environment, based on the [49][50][51] toolbox. Several segmentation metrics were calculated including accuracy (Equation (8)), sensitivity (Equation (9)), precision (Equation (10)), Matthews correlation coefficient (MCC) (Equation (11)), Dice index (Equation (12)), and Jaccard index (Equation (13)). The optimum results from this analysis were also confirmed with high-resolution ground truth data, provided after the digitization of an aerial orthophoto of the area, scale 1:5000.
where TP refers to true positive, TN to true negative, FN to false negative, and FP to false positive. The findings are discussed in the context of the archaeological site of "Nea Paphos". The overall methodology discussed earlier in this section is depicted in Figure 2.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 19 method. Then, the watershed algorithm is applied to the gradient image [46]. The watershed segmentation is a region-based method that has its origins in mathematical morphology [47]. In watershed segmentation, an image is regarded as a topographic landscape with ridges and valleys [48]. The segmentation performance, as a comparison of the segments produced from the 30 meter Landsat 8 image and the fused 10 meter Landsat 8 image was then examined. The first image was used as the ground truth image, while the second image was used as the slave image. This procedure was carried out once again in the MathWorks MATLAB R2016b environment, based on the [49][50][51] toolbox. Several segmentation metrics were calculated including accuracy (Equation (8)), sensitivity (Equation (9)), precision (Equation (10)), Matthews correlation coefficient (MCC) (Equation (11)), Dice index (Equation (12)), and Jaccard index (Equation (13)). The optimum results from this analysis were also confirmed with high-resolution ground truth data, provided after the digitization of an aerial orthophoto of the area, scale 1:5000.
where TP refers to true positive, TN to true negative, FN to false negative, and FP to false positive. The findings are discussed in the context of the archaeological site of "Nea Paphos". The overall methodology discussed earlier in this section is depicted in Figure 2.

Results
In this section, the fusion between the Landsat 8 (30 m) and the Sentinel 2 (10 m) images are demonstrated, while the results from the segmentation sensitivity follows.

Fusion Results
The general results of the fusion analysis between the Landsat 8 and Sentinel 2 images are shown in Figure 3. The first row of Figure 3, under the fused data label, presents the pansharpened Landsat 8 image based on the Gram-Schmidt method, the second row of Figure 3 presents the results after the application of the Brovey method, and the third and fourth rows indicate the pansharpened Landsat 8 image upon the application of the PCA and HSV methods, respectively. The first four columns of Figure 3, under the fused data label, present the different results of the above-mentioned pansharpening methods based on the four different Sentinel 2 (10 m) spectral bands used as the high-resolution input image. As mentioned earlier, these four spectral bands include the blue (Band 2), green (Band 3), red (Band 4), and NIR (Band 8) bands of Sentinel 2. The last column of Figure 3, under the raw data label, indicates the pseudo-color composites of the Landsat 8 (30 m) and Sentinel 2 (10 m) images at the red, green, and blue (RGB), and near-infrared red and green (NIR-R-G) visualization.   Compared to its original 30 m pixel resolution image, the fused pansharpened Landsat 8 image has an improved spatial resolution. The fusion is capable of capturing spatial details within the archaeological site of "Nea Paphos". This is true for all four different fusion methods applied in this study. It is worthwhile to notice that such details, within the archaeological site, were not visible in the original Landsat 30 m image (last column of Figure 3). The PCA fusion results are distinct as compared with the other fusion algorithms.
Moreover, while similar optical results are demonstrated for the visible bands (columns 1 to 3 of Figure 3, fused data) of each fusion method, the NIR band (fourth column of Figure 3) presents some differentiation. From visual interpretation, the Gram-Schmidt method seems to provide the most evident results as compared with the rest methods, which are similar with the original Sentinel 2 (10 m) image (shown in the last column of Figure 3). In this image, the excavated areas of the archaeological site (area c of Figure 1) are visible from visual interpretation, while the eastern boundary of the site with the modern city of Paphos, is becoming more distinct.
A closer look on the fusion results at specific areas of the archaeological site can be found in Figure 4 (for area b, Figure 1) and Figure 5 (for area c, Figure 1). These areas indicate the northern defensive wall of the city and other archaeological proxies, i.e., buried archaeological remains (area b, Figure 1), as well as existing excavated areas of the site (area c of Figure 1,).    Figure 1) using the four different pansharpening methods (Gram-Schmidt, Brovey, PCA, and HSV, at the first to fourth row, respectively) based on the different Sentinel 2 (10 m) spectral bands used as the panchromatic image (columns). The last column indicates the original 30 m Landsat 8 and Sentinel 2 10 m pixel resolution images at the red, green, and blue (RGB), and near-infrared green and blue (NIR-R-G) pseudo composites. The standing defensive wall is indicated with yellow arrows, and the archaeological proxies are indicated with white arrows.
As illustrated in Figure 4, the standing defensive wall (see yellow arrows at Figure 4), and archaeological proxies (see white arrows at Figure 4), are becoming visible in all outcomes of the fusion process as compared with the original Landsat 8 image. Similar observations, as those reported earlier, concerning the PCA fusion method are also valid. While it is difficult through visual interpretation to assess which method and spectral band better enhance the contrast of the archaeological proxies and standing defensive wall with the surrounding archaeological site, once again the Gram-Schmidt method seems to be the best fusion method.
Comparable observations with those of Figure 4, can also be reported for area c in Figure 5. The excavated areas of the archaeological site (see yellow arrows at Figure 5) as well as the eastern boundary of the site with the modern city of Paphos (see white arrows at Figure 5) are becoming visible after the application of the fusion methods, while the Gram-Schmidt method once again provides the clearest (sharpen) fused image.
Following the visual inspection of the fusion results (see Figures 3-5), a quantitative analysis evaluating the performance of each fusion method was carried out. The overall results of this analysis are grouped in Table 1. As mentioned in the previous section, five different quality image methods, namely the bias, image entropy, ERGAS, RASE, and RMSE, have been applied. Each method compares the spectral distance between the original (i.e., Landsat 8 image, 30 m) and the fused image (i.e., Landsat 8 fused image, 10 m), based on the different mathematical equations indicated in Equations (3) to (7). Results closest to zero indicate no significant spectral difference between the two images.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 19 row, respectively) based on the different Sentinel 2 (10 m) spectral bands used as the panchromatic image (columns). The last column indicates the original 30 m Landsat 8 and Sentinel 2 10 m pixel resolution images at the red, green, and blue (RGB), and near-infrared green and blue (NIR-R-G) pseudo composites. The standing defensive wall is indicated with yellow arrows, and the archaeological proxies are indicated with white arrows. Following the visual inspection of the fusion results (see Figures 3-5), a quantitative analysis evaluating the performance of each fusion method was carried out. The overall results of this analysis are grouped in Table 1. As mentioned in the previous section, five different quality image methods, namely the bias, image entropy, ERGAS, RASE, and RMSE, have been applied. Each method compares the spectral distance between the original (i.e., Landsat 8 image, 30 m) and the fused image (i.e., Landsat 8 fused image, 10 m), based on the different mathematical equations indicated in Equations (3) to (7). Results closest to zero indicate no significant spectral difference between the two images.
As presented in Table 1, the NIR spectral band for both Gram-Schmidt and Brovey transformations tends to give better results, with the exception of the Brovey method based on the entropy quality image method. In contrast, for the PCA and HSV fusion methods, the best results are obtained at the visible spectral bands, namely the blue, green, and red bands.
However, once we compare all the results, we can observe that the Gram-Schmidt method provides the best results, i.e., lower values for each different method, in relation to the rest of the fusion models. The score of the Gram-Schmidt method using the NIR spectral band of Sentinel 2 (Band 8) for the bias, entropy, ERGAS, RASE, and RMSE was estimated to 0.028, 0.360, 0.748, 3.271, As presented in Table 1, the NIR spectral band for both Gram-Schmidt and Brovey transformations tends to give better results, with the exception of the Brovey method based on the entropy quality image method. In contrast, for the PCA and HSV fusion methods, the best results are obtained at the visible spectral bands, namely the blue, green, and red bands.
However, once we compare all the results, we can observe that the Gram-Schmidt method provides the best results, i.e., lower values for each different method, in relation to the rest of the fusion models. The score of the Gram-Schmidt method using the NIR spectral band of Sentinel 2 (Band 8) for the bias, entropy, ERGAS, RASE, and RMSE was estimated to 0.028, 0.360, 0.748, 3.271, and 5.058, respectively, while the same score for the Brovey pansharpening method was calculated to 0.991, 1.654, 27.364, 100.783, and 236.473.
This finding is also compatible with the visual interpretation analysis and the outcomes of Figures 3-5, where the specific method seemed to provide the sharpest images for the archaeological site. In addition to this, and in alignment with the interpretation of the fused images at Figures 3-5, the best 10 m spectral band of the Sentinel 2 is the NIR band. Therefore, with this evaluation analysis over the archaeological site of "Nea Paphos", we conclude that the Gram-Schmidt fusion method based on the spectral Band 8 (NIR) of the Sentinel 2, better improves the visual interpretation of the 30 m Landsat 8 to that of a 10 m fused Landsat 8 image, without any significant spectral distortion of the original image.

Segmentation Analysis
According to these findings, we moved on to the next step of our analysis, investigating the sensitivity of the segmentation process of the fused image (Gram-Schmidt pansharpened Landsat 8 image) as compared with the original Landsat. The segmentation process was performed using different segmentation parameters, changing thus the "scale" parameter from 0 to 100, with a step of 10. The results from the segmentation process for area b of Figure 1 are shown in Figure 6 (for scales 10 to 50) and Figure 7 (for scales 60 to 100). Groups of pixels and segments are visualized in both images with a blue polygon.
Landsat 8 image is capable of capturing part of the defendsive wall and the archaeological proxies. The defensive wall is well captured at scales 60 and 70 (see Figure 7) since the segments capture various objects of interest at this area of the archaeological site, whereas for the next coarse scales (beyond the factor 70), the segments are quite extensive and not able to provide any useful information for this area (under-segmentation). Similar observations, not shown here, have also been reported for area c of Figure 1 (excavated archaeological areas).   Figure 7. The standing defensive wall is indicated with yellow arrows and the archaeological proxies are indicated with black arrows.
As presented in Figure 6, the original Landsat 8 for scales 10 to 40, do not capture the details of the archaeological site in this area, namely the standing defensive wall (indicated with yellow arrows in Figure 6) and the archaeological proxies (indicated with black arrows in Figure 6). The segmentation process at these scales provides almosts pixel-based segmentations (over-segmentation). Therefore the segments are too small and can not capture objects of interest.
Partially good segmentation of the original Landsat 8 is performed at a scale factor of 50 (last row of Figure 6), whereas groups pixels are aligned together, however, once again over-segmentation phenomenon can be also observed. Beyond this scale factor, i.e., for the scale factors 60 to 100 (see Figure 7), an under-segmentation phenomenon is visible since segments are becoming too large to map the objects of interest at the archaeological site.
Therefore, the segmentation process, with the original 30 m Landsat image tends to provide poor results, mainly due to its spatial resolution. In contrast, the fused Gram-Schmidt pansharpened Landsat 8 image is capable of capturing part of the defendsive wall and the archaeological proxies. The defensive wall is well captured at scales 60 and 70 (see Figure 7) since the segments capture various objects of interest at this area of the archaeological site, whereas for the next coarse scales (beyond the factor 70), the segments are quite extensive and not able to provide any useful information for this area (under-segmentation). Similar observations, not shown here, have also been reported for area c of Figure 1 (excavated archaeological areas).
Remote Sens. 2020, 12, x FOR PEER REVIEW  12 of 19 image (right). Groups of pixels (segments) are visualized in both images with a blue polygon. Similar results for scales 60 to 100 are shown in Figure 7. The standing defensive wall is indicated with yellow arrows and the archaeological proxies are indicated with black arrows.  Figure 6. The standing defensive wall is indicated with yellow arrows and the archaeological proxies are indicated with black arrows. Table 2 shows the results of the segmentation performance based on the seven different segmentation metrics, as described earlier in Equations (8) to (13). The segmentation performance of this Table is a Figure 1 using a scale factor from 60 to 100, with a step of 10, applied at the Gram-Schmidt pansharpened Landsat 8 image (left) and the original Landsat 8 image (right). Groups of pixels (segments) are visualized in both images with a blue polygon. Similar results for scales 10 to 50 are shown in Figure 6. The standing defensive wall is indicated with yellow arrows and the archaeological proxies are indicated with black arrows. Table 2 shows the results of the segmentation performance based on the seven different segmentation metrics, as described earlier in Equations (8) to (13). The segmentation performance of this Table is a result of the direct comparison of the Gram-Schmidt fused Landsat image and the  original Landsat image. As shown in the table, the scale factor 80 provides poor results for all metrics, whereas, for the scale factor 90, the segmentation of both images is the same (score equals to one). This is because the segmentation scale in both images is quite large, and therefore provides identical outcomes, leading an under-segmentation process. The most interesting results are found at scale factor 70, whereas the accuracy, sensitivity, and precision metrics are quite high (as also reported in the smaller scales) and the Dice and Jaccard indexes are the highest among all scales. Most important is that the Matthews correlation coefficient (MCC) score at the scale factor 70 (indicated by grey in Table 2) is the highest from all other scale factors. The MCC is an indicator of the quality of the segmentation as it takes into consideration all the scores from the true positives, true negatives, false positives, and false negatives results. When the value of MCC is one, this is an indicator of a perfect positive correlation between the segments. Conversely, when the classifier always misclassifies, then MCC gets a value of −1. Quantitative results presented in Table 2 are aligned with the visual interpretation and findings of Figures 6 and 7 reported earlier. These results suggest that the Gram-Schmidt fused Landsat image generates segments that are, on the one hand, aligned with the segments of the 30 m Landsat image (see Table 2), but, on the other hand, with more details, as evident in Figures 6 and 7. However, a comparison between the Gram-Schmidt fused Landsat image against ground truth data is needed. For this comparison, digitization of different areas of interest within the archaeological site of "Nea Paphos" was carried out in the ESRI ArcGIS v10.6 environment and used as a ground truth dataset. Digitized areas are shown as orange and yellow in Figure 8, based on the aerial orthophoto of the area (scale 1:5000) provided by the Department of Land and Surveys of Cyprus (see also Figure 1). These objects included both excavated large areas, visible in the center of Figure 8 with orange regular shapes, and other archaeological proxies indicated by yellow in Figure 8. The background of Figure 8 shows the segmentation of the Gram-Schmidt pansharpened Landsat 8 image using the NIR band of Sentinel. Using ground truth areas as input, it was found that scale factor 70 provided reasonably good results in regard to the segmentation based on both the Dice and Jaccard indices. The score for these two indices was estimated to be 0.687 and 0.523, respectively, which was considered reasonable due to the spatial resolution (10 m) of the fused image as compared with the very high detail (~1 m error) of the orthophoto.
Important details within the archaeological site of "Nea Paphos" are indicated with yellow arrows and areas a-f, in both The segmentation difference between the original Landsat 8 image and the fused Gram-Schmidt pansharpened image, for a scale factor 70, is shown in Figure 10.  The segmentation difference between the original Landsat 8 image and the fused Gram-Schmidt pansharpened image, for a scale factor 70, is shown in Figure 10.

Conclusions
The segmentation process is of great importance to the object-oriented analysis of multispectral satellite images. This paper aims to investigate and evaluate the impact of the fusion of freely distributed satellite-based images such as those of Landsat 8 and Sentinel 2, to improve, at a first step, the spatial resolution of the first image and, consequently, to evaluate the segmentation performance of the fused images based on different scale factors.
The Gram-Schmidt pansharpening method, the near-infrared spectral 10 meter band of Sentinel 2 (namely Band 8), and the scale factors at approximately 50 to 70 were the best indicators for the segmentation process for the archaeological site of "Nea Paphos". This conclusion was supported both from the visual interpretation of the overall results, as well as from the quantitative analysis carried out and demonstrated in Tables 1 and 2. While a single segmentation method was implemented, here, other segmentation algorithms and strategies could be applied to improve the overall performance of the fused datasets. In our study, the watershed transformation and Sobel filter were used for the segmentation process in various scale factors, which are, however, strongly sensitive to image noise. Therefore, future research could focus in this direction, and therefore  Figure 9 shows the segmentation of the Gram-Schmidt method at the scale factor 70 over the high-resolution aerial orthophoto image (Figure 9, left) and over the segmented Landsat 30 m pixel resolution (Figure 9, right). As demonstrated, the segmentation at this range effectively captures the defensive wall at the northern part of the archaeological site (arrow a), as well as other areas within the "Nea Paphos" site. The fusion of the Landsat 8 image improves the segmentation process as compared with the direct segmentation of the 30 m Landsat 8 image at the same scale factor. Large segments observed in the Landsat 8 (Figure 9b), which removes details and geometrical shapes of the archaeological site, can now be partially recovered in the Gram-Schmidt fused Landsat image. The segmentation difference between the original Landsat 8 image and the fused Gram-Schmidt pansharpened image, for a scale factor 70, is shown in Figure 10.

Conclusions
The segmentation process is of great importance to the object-oriented analysis of multispectral satellite images. This paper aims to investigate and evaluate the impact of the fusion of freely distributed satellite-based images such as those of Landsat 8 and Sentinel 2, to improve, at a first step, the spatial resolution of the first image and, consequently, to evaluate the segmentation performance of the fused images based on different scale factors.
The Gram-Schmidt pansharpening method, the near-infrared spectral 10 m band of Sentinel 2 (namely Band 8), and the scale factors at approximately 50 to 70 were the best indicators for the segmentation process for the archaeological site of "Nea Paphos". This conclusion was supported both from the visual interpretation of the overall results, as well as from the quantitative analysis carried out and demonstrated in Tables 1 and 2. While a single segmentation method was implemented, here, other segmentation algorithms and strategies could be applied to improve the overall performance of the fused datasets. In our study, the watershed transformation and Sobel filter were used for the segmentation process in various scale factors, which are, however, strongly sensitive to image noise. Therefore, future research could focus in this direction, and therefore minimize the impact of noise (e.g., due to radiometric or atmospheric effects) and improve the segmentation analysis.
Although the results reported, in this study, are only representative of the specific archaeological site of interest, the overall methodological framework could be adopted in any other archaeological site where Landsat and Sentinel images can be found. Important aspects of these datasets are their open access policy distribution, with no acquisition cost, and high temporal revisit times (16 days for Landsat 8 and five days for Sentinel 2 sensors) that are ideal for capturing archaeological sites and landscapes in specific time windows (upon cloud coverage). The proposed methodology is only applicable, however, for images acquired after 2015, when the Sentinel optical sensors datasets (2A and 2B) are available along with the Landsat 8 images. Since Sentinel 2 Level 1C datasets are not atmospherically corrected, this needs to be taken into consideration prior to the fusion. Images taken with minimum cloud coverage, on non-hazy days, are expected to have minimum atmospheric impact, and this can be evaluated using simple image-based atmospheric correction methods such as those of darkest pixel.
Segmentation of Landsat fused images using Sentinel (10 m) bands are expected to better support object-oriented (OBIA) classifications and increase the detection rate of archaeological proxies, therefore, expanding the scope of raw 30 m resolution Landsat images.
Future work is expected to focus in this direction, aiming to further study synergistic use of different sensors, integrated to extract better results for satellite-based archaeological prospection and management, as well as work in other periods of times, whereas archaeological proxies can be better enhanced [52].