Working with Gaussian Random Noise for Multi-Sensor Archaeological Prospection: Fusion of Ground Penetrating Radar Depth Slices and Ground Spectral Signatures from 0.00 m to 0.60 m below Ground Surface

: The integration of di ﬀ erent remote sensing datasets acquired from optical and radar sensors can improve the overall performance and detection rate for mapping sub-surface archaeological remains. However, data fusion remains a challenge for archaeological prospection studies, since remotely sensed sensors have di ﬀ erent instrument principles, operating in di ﬀ erent wavelengths. Recent studies have demonstrated that some fusion modelling can be achieved under ideal measurement conditions (e.g., simultaneously measurements in no hazy days) using advance regression models, like those of the nonlinear Bayesian Neural Networks. This paper aims to go a step further and investigate the impact of noise in regression models, between datasets obtained from ground-penetrating radar (GPR) and portable ﬁeld spectroradiometers. Initially, the GPR measurements provided three depth slices of 20 cm thickness, starting from 0.00 m up to 0.60 m below the ground surface while ground spectral signatures acquired from the spectroradiometer were processed to calculate 13 multispectral and 53 hyperspectral indices. Then, various levels of Gaussian random noise ranging from 0.1 to 0.5 of a normal distribution, with mean 0 and variance 1, were added at both GPR and spectral signatures datasets. Afterward, Bayesian Neural Network regression ﬁtting was applied between the radar (GPR) versus the optical (spectral signatures) datasets. Di ﬀ erent regression model strategies were implemented and presented in the paper. The overall results show that fusion with a noise level of up to 0.2 of the normal distribution does not dramatically drop the regression model between the radar and optical datasets (compared to the non-noisy data). Finally, anomalies appearing as strong reﬂectors in the GPR measurements, continue to provide an obvious contrast even with noisy regression modelling.

The paper presents the overall methodology adopted as well as the results obtained from the various simulations working with noisy data. Initially, a short description of the case study area and the datasets are presented, following by the methodology section. Then the overall results of this simulation study are provided, following by a general discussion on the key findings of the study.

Case Study Area
One of the Hungarian National Parks, namely the Vésztő-Mágor Tell, situated in the southeastern Great Hungarian Plain, was selected as the case study ( Figure 1). The site, which approximately covers 4.25 hectares at a height of 9 m a.s.l., was investigated in the previous years as part of the Körös Regional Archaeological Project [28] with various non-invasive techniques including magnetics, GPR, electrical resistivity tomography (ERT) and Electromagnetic induction (EMI) methods. Previous excavations during the 1970s and 1980s revealed that the area was initially settled by Szakálhát culture farmers during the late Middle Neolithic period, while the area was systematically occupied until the Late Neolithic period (c. 5000-4600 B.C. calibrated). On top of this stratigraphy, the layers found suggests the abandonment of the site until the formation of an Early Copper Age settlement and later on by a Middle Copper Age settlement. During the 11th century AD, a church followed by a monastery was established on top of the Tell. After various modifications and reconstructions, the monastery was destroyed during the early 19th century. Currently, in addition to the archaeological museum located within the historic wine cellar, the 1986 excavation trench can be seen with several features from different periods remaining in situ. More information regarding the archaeological context of the Vésztő-Mágor Tell can be found in [29][30][31][32][33].

Data Sources
In the last phases of the Körös Regional Archaeological Project (KRAP), which was directed by William A. Parkinson of the Department of Anthropology, The Field Museum of Chicago, and Attila Gyucha of the Field Service of Cultural Heritage, Hungary, two different ground prospection methods were implemented in addition to the previous ground-based prospection techniques. The first one was a traditional GPR survey using a Noggin Plus (Sensors & Software, Ontario, Canada) GPR with a 250 MHz antenna while the second method was based on a portable field spectroradiometer, namely the GER-1500 spectroradiometer (Spectra Vista Corporation, New York, NY, USA), with a spectral range of 350-1050 nm. The area of measurements is shown in Figure 1 above. The red polygon indicates the area covered from the ground spectroradiometric measurements while the blue polygon indicates the areas covered from the GPR geophysical survey.

Methodology
The overall methodology of this study is presented in Figure 2. Both spectroradiometric and GPR measurements were collected moving along parallel transects 0.5 m apart, with 5 cm sampling along the lines for the GPR, covering an area of 20 × 90 m (1800 square meters). The data have been collected simultaneously to eliminate the influence of noise from other external factors (mainly climatic changes and variation of soil conditions, such as humidity). Once the data collection was accomplished, pre-processing steps including noise removal and filtering were applied. Then depth slices of 20 cm thickness were produced from the GPR datasets. The ground spectroradiometric measurements were resampled to multispectral bands (blue; green; red; near-infrared) using a relative spectral response (RSR) filter (see Equation (1)). The RSR filters define the instrument's relative sensitivity to radiance in various parts of the electromagnetic spectrum. RSR filters have a value range from 0 to 1 and are unit-less since they are relative to the peak response.
where R band = reflectance at a range of wavelength (e.g., blue band); Ri = reflectance at a specific wavelength (e.g., R 450 nm) and RSRi = relative response value at the specific wavelength. Based on the multispectral bands and the hyperspectral bands of the spectroradiometer, more than 65 hyperspectral and multispectral indices were calculated. The overall equations and references of the vegetation indices can be found in Appendix A. From previous studies [25,26], the first three depth slices (i.e., 0.00-0.20; 0.20-0.40 and 0.40-0.60 m below ground surface) were correlated with the optical products taken from the spectroradiometer (i.e., vegetation indices). Depths beyond the 0.60 m below ground surface did not provide any reliable results [25]. More details regarding data processing can be found in [25,26].
In our example, a randomly normally distributed noise of a standard deviation of 0.1; 0.2; 0.3; 0.4 and 0.5 (hereafter referred as the noise of level 0.1; 0.2; 0.3; 0.4 and 0.5 respectively) was gradually added at the three GPR depth slices and the vegetation indices. Various assessments were carried out based on all levels of noise aiming to understand the correlation performance of the radar datasets against the optical products. This correlation, as mentioned earlier, was based on the Bayesian Neural Network regression fitting, following the same procedure as in Agapiou and Sarris [25].
The neural network fitting was processed by the MathWorks MATLAB R2016b software. The Bayesian Neural Network is a neural network with a prior distribution on its weights, where the likelihood for each data is given by Equation (2). p y n w, x n , σ 2 = Normal y n NN(; w), σ 2 (2) where NN is a neural network whose weights and biases originate from the latent variable w [34]. From the whole dataset, 85% was used for training purposes while the remaining 15% was used for testing the fusion models. The results from the above regression were compared against the real GPR datasets and the non-noisy regression models.

Figure 2.
Overall methodology followed in this study.

Results
The various results produced during this study are presented below. The results are grouped into four categories for better interpretation. Section 4.1 presents the impact of the random noise added at the multispectral bands regressed against the GPR depth slices. The next section (Section 4.2) includes the results from the regression models between the multispectral bands against the noisy GPR depth slices (from 0.00 to 0.60 m), while the impact of the random noise added at both multispectral bands and the GPR depth slices are presented in Section 4.3. Finally, the impact of the random noise added to vegetation indices against the GPR depth slices is demonstrated in Section 4.4.

Random Noise at Spectral Bands
Normally distributed random noise with a standard deviation of 0.1; 0.2; 0.3; 0.4 and 0.5 was added in the four multispectral bands (Band 1: Blue; Band 2: Green; Band 3: Red and Band 4: Nearinfrared) after the use of the RSR filters as shown in Figure 3. Column (a) shows the initial spectral bands as extracted from the spectroradiometric campaigns following a kriging interpolation method. Column (a) is considered in this case as the ground truth dataset for comparison purposes. Columns (b) to (f) shows the multispectral bands after the inclusion of the random noise starting from a level of 0.1 up to 0.5, respectively. The first row of Figure 3 shows the results for the first band (blue), while the second (green), third (red) and the fourth (near-infrared) bands are demonstrated on the second, third and fourth rows of Figure 3, respectively. It is obvious that the noise added in the spectral bands alters the final interpretation compared to the truth datasets (Column (a)), especially in the nearinfrared part of the spectrum (see 4th row). However, strong anomalies-as depicted from the spectroradiometer-as those observed at around 80 m on the north axis (see the arrow of Figure 3) are still visible even in datasets with a high level of noise.
Preliminary coring results suggest that such linear anomalies, like those observed in Figure 3 (see black arrow) should be linked with ditches, while some of them (not shown here) may exceed 4 m in depth. Structural similarities with other Neolithic fortifications led archaeologists to suspect that part of the ditch and palisade system at the Vésztő-Mágor Tell was established during this period [28].

Results
The various results produced during this study are presented below. The results are grouped into four categories for better interpretation. Section 4.1 presents the impact of the random noise added at the multispectral bands regressed against the GPR depth slices. The next section (Section 4.2) includes the results from the regression models between the multispectral bands against the noisy GPR depth slices (from 0.00 to 0.60 m), while the impact of the random noise added at both multispectral bands and the GPR depth slices are presented in Section 4.3. Finally, the impact of the random noise added to vegetation indices against the GPR depth slices is demonstrated in Section 4.4.

Random Noise at Spectral Bands
Normally distributed random noise with a standard deviation of 0.1; 0.2; 0.3; 0.4 and 0.5 was added in the four multispectral bands (Band 1: Blue; Band 2: Green; Band 3: Red and Band 4: Near-infrared) after the use of the RSR filters as shown in Figure 3. Column (a) shows the initial spectral bands as extracted from the spectroradiometric campaigns following a kriging interpolation method. Column (a) is considered in this case as the ground truth dataset for comparison purposes. Columns (b) to (f) shows the multispectral bands after the inclusion of the random noise starting from a level of 0.1 up to 0.5, respectively. The first row of Figure 3 shows the results for the first band (blue), while the second (green), third (red) and the fourth (near-infrared) bands are demonstrated on the second, third and fourth rows of Figure 3, respectively. It is obvious that the noise added in the spectral bands alters the final interpretation compared to the truth datasets (Column (a)), especially in the near-infrared part of the spectrum (see 4th row). However, strong anomalies-as depicted from the spectroradiometer-as those observed at around 80 m on the north axis (see the arrow of Figure 3) are still visible even in datasets with a high level of noise.   Preliminary coring results suggest that such linear anomalies, like those observed in Figure 3 (see black arrow) should be linked with ditches, while some of them (not shown here) may exceed 4 m in depth. Structural similarities with other Neolithic fortifications led archaeologists to suspect that part of the ditch and palisade system at the Vésztő-Mágor Tell was established during this period [28].
In order to quantify the impact of the added noise to the regression model, all four multispectral bands (in all levels of noise, including level 0-noiseless) were regressed separately against the three GPR depth slices. The Pearson regression value (r), as well as the relative Mean Square Error (MSE) compared with level 0 (noiseless datasets), used in the next sections, are criteria of the regression performance and are presented in Equations (3) and (4), respectively: where x is the observed GPR-measurement, y is the BNN modelled GPR result, and n is the total number of measurements. Figure 4a shows the results of the Pearson regression value and Figure 4b the results from the MSE analysis. We can observe that the added noise decreases steadily the correlation value r while in contrary the MSE increases progressively for all depth slices. The correlation coefficient r falls from over 0. In order to quantify the impact of the added noise to the regression model, all four multispectral bands (in all levels of noise, including level 0-noiseless) were regressed separately against the three GPR depth slices. The Pearson regression value (r), as well as the relative Mean Square Error (MSE) compared with level 0 (noiseless datasets), used in the next sections, are criteria of the regression performance and are presented in Equations (3) and (4), respectively: where x is the observed GPR-measurement, y is the BNN modelled GPR result, and n is the total number of measurements. Figure

Random Noise at GPR Depth Slices
In a similar way as before, the Bayesian Neural Network correlation was applied after adding random noise only to the GPR depth slices. Column (a) of Figure 5, shows the initial raw GPR data (noiseless datasets), while the rest of the columns (b to f) shows the simulated GPR depth slices with a noise level from 0.1 to 0.5, respectively. Results from the depth slice 0.00-0.20 m are demonstrated in the first row, while the depth slices 0.20-0.40 and 0.40-0.60 m are shown in the second and thirdrow correspondingly. From the visual interpretation of the results, it becomes evident that noise of levels 0.1 and 0.2 share similar patterns of anomalies as the one depicted at the original datasets (see first column of Figure 5, noiseless GPR depth slices), while noise beyond the level of 0.2 saturates the overall GPR signal. However, once again, the strong linear anomaly at around 80 m on the y-axis remains visible in all levels of noise (see arrow in Figure 5). Figure 6 shows the statistics from this Bayesian Neural Network regression fitting of the multispectral bands against the GPR depth slices (with and without noise). The results are shown once again in relation to the Pearson correlation coefficient (r) and the relative MSE in terms of percentage compared with the noiseless datasets. The correlation analysis shows a comparable trend as observed in the previous experiment (see Section 4.1; Figure 4): the correlation coefficient value drops gradually with the increase of noise for all depth slices, having similar linear trend, while the MSE increases exponentially with the increase of noise level. It should be mentioned that for noise levels 0.1 and 0.2 the correlation coefficient decreases slightly from 0.50 to 0.40 for the first depth slice, from 0.42 to 0.38 for the second depth slice and from 0.48 to 0.46 for the third depth slice. Even though MSE for these depth slices increases, it is worth to notice the large difference in the MSE between the first depth slice (0.00-0.20m) and the rest of the depth slices.

Random Noise at GPR Depth Slices
In a similar way as before, the Bayesian Neural Network correlation was applied after adding random noise only to the GPR depth slices. Column (a) of Figure 5, shows the initial raw GPR data (noiseless datasets), while the rest of the columns (b to f) shows the simulated GPR depth slices with a noise level from 0.1 to 0.5, respectively. Results from the depth slice 0.00-0.20 m are demonstrated in the first row, while the depth slices 0.20-0.40 and 0.40-0.60 m are shown in the second and third-row correspondingly. From the visual interpretation of the results, it becomes evident that noise of levels 0.1 and 0.2 share similar patterns of anomalies as the one depicted at the original datasets (see first column of Figure 5, noiseless GPR depth slices), while noise beyond the level of 0.2 saturates the overall GPR signal. However, once again, the strong linear anomaly at around 80 m on the y-axis remains visible in all levels of noise (see arrow in Figure 5). Figure 6 shows the statistics from this Bayesian Neural Network regression fitting of the multispectral bands against the GPR depth slices (with and without noise). The results are shown once again in relation to the Pearson correlation coefficient (r) and the relative MSE in terms of percentage compared with the noiseless datasets. The correlation analysis shows a comparable trend as observed in the previous experiment (see Section 4.1; Figure 4): the correlation coefficient value drops gradually with the increase of noise for all depth slices, having similar linear trend, while the MSE increases exponentially with the increase of noise level. It should be mentioned that for noise levels 0.1 and 0.2 the correlation coefficient decreases slightly from 0.50 to 0.40 for the first depth slice, from 0.42 to 0.38 for the second depth slice and from 0.48 to 0.46 for the third depth slice. Even though MSE for these depth slices increases, it is worth to notice the large difference in the MSE between the first depth slice (0.00-0.20m) and the rest of the depth slices.

Random Noise at Spectral Bands and GPR Depth Slices
Afterward, random noise was added at both multispectral bands and the GPR depth slices starting from level 0.1 and increasing by 0.1 until the level of 0.5. Again, the (noisy) spectral bands were set as inputs and the (noisy) GPR depth slices were set as targets in the Bayesian Neural Network regression fitting models. Once the correlation was estimated then the optical data (i.e., spectral bands) were simulated as GPR depth slices as shown in Figure 7.

Random Noise at Spectral Bands and GPR Depth Slices
Afterward, random noise was added at both multispectral bands and the GPR depth slices starting from level 0.1 and increasing by 0.1 until the level of 0.5. Again, the (noisy) spectral bands were set as inputs and the (noisy) GPR depth slices were set as targets in the Bayesian Neural Network regression fitting models. Once the correlation was estimated then the optical data (i.e., spectral bands) were simulated as GPR depth slices as shown in Figure 7.   Figure 7 is that the linear anomaly that is depicted in all depth slices (see arrow in Figure 7) continues to be visible. The high contrast of this linear anomaly is noticeable in all scenarios, even with a noise level of 0.2, allowing thus to suggest that some correlation could be achieved for strong anomalies, even in cases where no ideal conditions of measurements are feasible.
The overall statistics from the above-mentioned analysis are shown in Figure 8. The correlation coefficient value r ranges for noise level 0.1 from 0.30 up to 0.45 for all depth slices. However, the correlation coefficient value r decreases even further for the noisiest datasets (level of 0.2 at spectral bands and level 0.1 at the GPR depth slices). The correlation achieved by Agapiou and Sarris [25] is shown as with a dot in Figure 8 for each depth slice. The difference achieved in the correlation using the truth data (indicated with dots) and the noisy datasets of level 0.2 is obvious in Figure 8.
Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 23 0.1 noise at spectral band and 0.1 noise at raw GPR data; c-I: Simulated GPR data at depth slice of 0.00-0.20 m below ground surface with 0.2 noise at spectral band and 0.1 noise at raw GPR data; c-II: Simulated GPR data at depth slice of 0.20-0.40 m below ground surface with 0.2 noise at spectral band and 0.1 noise at raw GPR data and c-III: Simulated GPR data at depth slice of 0.40-0.60 m below ground surface with 0.2 noise at spectral band and 0.1 noise at raw GPR data.
Column (a) of Figure 7 shows the raw GPR data at the depth slice of 0.00-0.20; 0.20-0.40 and 0.40-0.60 m below ground surface. Column (b) of Figure 7, shows the simulated GPR data (from the optical data) at depth slice of 0.00-0.20, 0.20-0.40 and 0.40-0.60 m below ground surface with the noise of level 0.1. The last column (c) of Figure 7, shows the simulated results as before, with a noise level of 0.2. The most interesting result from Figure 7 is that the linear anomaly that is depicted in all depth slices (see arrow in Figure 7) continues to be visible. The high contrast of this linear anomaly is noticeable in all scenarios, even with a noise level of 0.2, allowing thus to suggest that some correlation could be achieved for strong anomalies, even in cases where no ideal conditions of measurements are feasible.
The overall statistics from the above-mentioned analysis are shown in Figure 8. The correlation coefficient value r ranges for noise level 0.1 from 0.30 up to 0.45 for all depth slices. However, the correlation coefficient value r decreases even further for the noisiest datasets (level of 0.2 at spectral bands and level 0.1 at the GPR depth slices). The correlation achieved by Agapiou and Sarris [25] is shown as with a dot in Figure 8 for each depth slice. The difference achieved in the correlation using the truth data (indicated with dots) and the noisy datasets of level 0.2 is obvious in Figure 8.

Random Noise and Vegetation Indices
Even though a moderate correlation coefficient r was achieved in the previous experiments, these results are aligned with the outcomes presented by Agapiou and Sarris [15]. The correlation is expected to increase once we use vegetation indices (both multispectral and hyperspectral indices) in the regression model. For this reason, the regression analysis was extended by analyzing the correlation performance against the (a) NDVI, (b) 13 multispectral and (c) 53 hyperspectral indices, listed in Tables A1 and A2 (Appendix A).
The results of the regression modelling between the NDVI, the multispectral and hyperspectral indices against all the noisy GPR depth slices of level 0.1 and level 0.2 are shown in Figure 9a

Random Noise and Vegetation Indices
Even though a moderate correlation coefficient r was achieved in the previous experiments, these results are aligned with the outcomes presented by Agapiou and Sarris [15]. The correlation is expected to increase once we use vegetation indices (both multispectral and hyperspectral indices) in the regression model. For this reason, the regression analysis was extended by analyzing the correlation performance against the (a) NDVI, (b) 13 multispectral and (c) 53 hyperspectral indices, listed in Tables A1 and A2 (Appendix A).
The results of the regression modelling between the NDVI, the multispectral and hyperspectral indices against all the noisy GPR depth slices of level 0.1 and level 0.2 are shown in Figure 9a,b, respectively. The correlation coefficient score is much higher-up to 0.61-compared to the previous outcomes discussed in . As in earlier studies (see [25,26]), the hyperspectral indices tend to provide the highest correlation coefficient values ranging from 0.50-0.60 for all depth slices. The regression score slightly decreases in the noisiest datasets (Figure 9b) for all combinations of correlation.
respectively. The correlation coefficient score is much higher-up to 0.61-compared to the previous outcomes discussed in sections 4.1-4.3. As in earlier studies (see [25,26]), the hyperspectral indices tend to provide the highest correlation coefficient values ranging from 0.50-0.60 for all depth slices. The regression score slightly decreases in the noisiest datasets (Figure 9b Tables A1 and A2 respectivelly, Appendix A). Correlation achieved with the noiseless dataset is shown with a dot for each depth slice.
A similar experiment was also carried out with noisy multispectral bands of levels 0.1 and 0.2, which subsequently were used to estimate 13 multispectral indices (see Table A1, Appendix A). The added noise dropped the correlation coefficient value from 0.50 to 0.40 and then to 0.35 for the first depth slice, 0.42 to 0.38 and then to 0.18 for the second depth slice and 0.50 to 0.44 and then to 0.38 for the last depth slice, as estimate for noiseless, level 0.1 and level 0.2 noise, respectively. The overall A similar experiment was also carried out with noisy multispectral bands of levels 0.1 and 0.2, which subsequently were used to estimate 13 multispectral indices (see Table A1, Appendix A). The added noise dropped the correlation coefficient value from 0.50 to 0.40 and then to 0.35 for the first depth slice, 0.42 to 0.38 and then to 0.18 for the second depth slice and 0.50 to 0.44 and then to 0.38 for the last depth slice, as estimate for noiseless, level 0.1 and level 0.2 noise, respectively. The overall results are quite similar to the correlation achieved by Agapiou and Sarris [25] and illustrated as dots in Figure 10. This allows us to assume that the impact of noise up to a level of 0.2 does not dramatically change the correlation fitting score (r) between the optical and the radar datasets.
results are quite similar to the correlation achieved by Agapiou and Sarris [25] and illustrated as dots in Figure 10. This allows us to assume that the impact of noise up to a level of 0.2 does not dramatically change the correlation fitting score (r) between the optical and the radar datasets.  Table A1, Appendix A). Correlation achieved with the noiseless dataset is shown with a dot for each depth slice.

Discussion
As was expected, the added noise to the datasets changed the overall regression fitness of the models. The previous section provided four different examples of how random noise may influence the overall regression between the radar and the optical ground datasets. The overall findings suggest that fusion between diverse datasets, as those described here, can be challenging and further research is needed to improve the overall performance of the models by using more sophisticated training models. Based on the results presented before, we can conclude the following four general findings which are summarized below.


Finding-1: Noise applied either in the optical or the radar datasets decrease the overall correlation fitting performance. However, the noise of levels 0.1 and 0.2 does not dramatically drop the regression value, especially in the case when we used hyperspectral indices in the regression models. This noise can be for instance generated from atmospheric effects or from the sensitivity level of the sensor. Therefore, any potential noise can be minimized through a series of pre-processing steps as, for instance, the radiometric and atmospheric corrections implemented to optical datasets. In addition, external noise can be minimized once the datasets are collected concurrently, or at least when the climate conditions are the same (e.g., no humidity; no rainfall in the previous days). A survey protocol that minimizes potential noise can be of great importance in these cases as it will provide the overall guidelines and framework for the data collection in the field.  Finding-2: Anomalies appearing as strong reflectors in the GPR/spectroradiometer measurements, continue to provide an obvious contrast even with noisy datasets. Therefore, the synergistic use of diverse datasets can be implemented for detecting strong anomalies, which are usually identified as priority areas of interest in archaeological surveys. Strong anomalies are expected to be found in areas where the sub-surface target (archaeological remains?) has different properties in contrast to the surrounding soil matrix. The detection or not of such anomalies is also based on the capabilities and characteristics of the remote sensors used in the survey, which are sensitive to only a part of the spectrum.  Table A1, Appendix A). Correlation achieved with the noiseless dataset is shown with a dot for each depth slice.

Discussion
As was expected, the added noise to the datasets changed the overall regression fitness of the models. The previous section provided four different examples of how random noise may influence the overall regression between the radar and the optical ground datasets. The overall findings suggest that fusion between diverse datasets, as those described here, can be challenging and further research is needed to improve the overall performance of the models by using more sophisticated training models. Based on the results presented before, we can conclude the following four general findings which are summarized below.

•
Finding-1: Noise applied either in the optical or the radar datasets decrease the overall correlation fitting performance. However, the noise of levels 0.1 and 0.2 does not dramatically drop the regression value, especially in the case when we used hyperspectral indices in the regression models. This noise can be for instance generated from atmospheric effects or from the sensitivity level of the sensor. Therefore, any potential noise can be minimized through a series of pre-processing steps as, for instance, the radiometric and atmospheric corrections implemented to optical datasets. In addition, external noise can be minimized once the datasets are collected concurrently, or at least when the climate conditions are the same (e.g., no humidity; no rainfall in the previous days). A survey protocol that minimizes potential noise can be of great importance in these cases as it will provide the overall guidelines and framework for the data collection in the field. • Finding-2: Anomalies appearing as strong reflectors in the GPR/spectroradiometer measurements, continue to provide an obvious contrast even with noisy datasets. Therefore, the synergistic use of diverse datasets can be implemented for detecting strong anomalies, which are usually identified as priority areas of interest in archaeological surveys. Strong anomalies are expected to be found in areas where the sub-surface target (archaeological remains?) has different properties in contrast to the surrounding soil matrix. The detection or not of such anomalies is also based on the capabilities and characteristics of the remote sensors used in the survey, which are sensitive to only a part of the spectrum. • Finding-3: Vegetation indices provided the best correlation between the optical and GPR datasets. Figure 11 summarizes the overall results regarding the regression values for all different regression combinations between the spectroradiometric and GPR datasets, with an added noise of levels 0.1 and 0.2. The multispectral and hyperspectral indices tend to give the highest regression values r of up to 0.6 while the remaining datasets provided lower regression values (from 0.1 up to 0.5). This observation was compatible with the results found by Agapiou and Sarris in previous studies [25], indicating that processing of the raw datasets (e.g., reflectance values into vegetation indices) can further enhance the contrast of the sub-surface target. While the use of multispectral indices can be generated by most of the earth observation sensors, the use of hyperspectral indices is restricted to a limited number of sensors such as EO-Hyperion 1 (not active), and the Environmental Mapping and Analysis Program (EnMAP). • Finding-4: Regression fitting beyond the 0.6 meters (not shown in this paper) provided very low correlation coefficient (less than 0.10), allowing us to suggest that any correlation between optical and GPR can be performed only for the upper layers of soil until a depth of approximately 0.6 m. This was once again compatible with the earlier studies of the authors [25]. However, this finding needs to be further investigated with other types of vegetation having different root system characteristics and under different soil properties.
Remote Sens. 2019, 11, x FOR PEER REVIEW 15 of 23  Finding-3: Vegetation indices provided the best correlation between the optical and GPR datasets. Figure 11 summarizes the overall results regarding the regression values for all different regression combinations between the spectroradiometric and GPR datasets, with an added noise of levels 0.1 and 0.2. The multispectral and hyperspectral indices tend to give the highest regression values r of up to 0.6 while the remaining datasets provided lower regression values (from 0.1 up to 0.5). This observation was compatible with the results found by Agapiou and Sarris in previous studies [25], indicating that processing of the raw datasets (e.g., reflectance values into vegetation indices) can further enhance the contrast of the sub-surface target. While the use of multispectral indices can be generated by most of the earth observation sensors, the use of hyperspectral indices is restricted to a limited number of sensors such as EO-Hyperion 1 (not active), and the Environmental Mapping and Analysis Program (EnMAP).  Finding-4: Regression fitting beyond the 0.6 meters (not shown in this paper) provided very low correlation coefficient (less than 0.10), allowing us to suggest that any correlation between optical and GPR can be performed only for the upper layers of soil until a depth of approximately 0.6 m. This was once again compatible with the earlier studies of the authors [25]. However, this finding needs to be further investigated with other types of vegetation having different root system characteristics and under different soil properties. Even though that the previous findings are solely based on the experiments carried out in the archaeological site of the Hungarian Plain, and the specific measurements carried out in that ground survey (GPR and portable field spectroradiometer), the findings can be projected to other case studies that share similar climatic and soil sub-surface properties. This is something that needs, however, to be further examined and studied not only in different archaeological sites but also using a different type of sensors. While fluctuations of regression fitting models can be recorded in the future working in other archaeological landscapes or working with different sensors, the overall findings could be used so as to develop new strategies of how we can combine diverse archaeological prospection datasets, moving beyond the traditional layer super-imposed methods. Even though that the previous findings are solely based on the experiments carried out in the archaeological site of the Hungarian Plain, and the specific measurements carried out in that ground survey (GPR and portable field spectroradiometer), the findings can be projected to other case studies that share similar climatic and soil sub-surface properties. This is something that needs, however, to be further examined and studied not only in different archaeological sites but also using a different type of sensors. While fluctuations of regression fitting models can be recorded in the future working in other archaeological landscapes or working with different sensors, the overall findings could be used so as to develop new strategies of how we can combine diverse archaeological prospection datasets, moving beyond the traditional layer super-imposed methods.

Conclusions
Remote sensing technologies are nowadays well adopted for archaeological prospection studies focused on the detection of concealed archaeological remains beneath the ground surface. However, as stated in [35], "there is no one all-purpose remote sensing dataset on which the archaeologist can rely that will uncover all evidence of human occupation". Therefore, the need for synergistic use of different remote sensors, from space, air, and ground, is essential for the better understanding of the archaeological context.
The exploitation of diverse multi-sensory technology can improve the overall detection performance of sub-surface targets with respect to recall and precision. While in the last few years, several studies have been conducted using various sensors, either optical or radar, research regarding the fusion of these datasets is still very limited in the archaeological prospection domain.
As recently stated by Luo et al. [35], fusion of ground and spaceborne remote sensing data is one of the main barriers to the realization of integrated space-ground applications. To overcome this barrier, the different types of remote sensing datasets need to be harmonized and fused together. Once this achieved, then new ways of archaeological prospections are foreseen, linking ground source datasets with satellite observations, bringing together the advantages of each method for the benefit of the archaeological prospection.
The need to develop novel fusion techniques of different sensors is expected to be increased in the future due to the ongoing developments of multi-sensor equipment. These technological improvements can give higher spatial accuracy datasets, with even new wavelengths of observations (as the case of satellite radar sensors).
In the recent past, research studies [25,26] have shown that some correlation between optical and radar datasets for the upper layers of soil (up to 0.60 m beneath surface) can be achieved if we explore advance fitting regression models like those provided by the Bayesian neural networks. In addition, the regression fitness is improved once we process the initial datasets as the case of the reflectance values to vegetation indices.
A critical aspect of the previous studies, i.e., [25,26] was the fact that the datasets used between the ground spectroradiometer and the GPR measurements were taken simultaneously; thus, minimising potential noise. So, in any other case, the regression model would suffer from addition noise due to the impact on the environment. However, as shown in this study, even in the case where ideal conditions (i.e., simultaneously measurements) are not met, the impact of noise (in terms of standard deviation levels of 0.1 and 0.2 of a Gaussian distribution) can still provide reliable results. Still, as expected, the regression fitness is lower than the regression performed in ideal conditions. Future studies will explore the benefits from these outcomes such as the correlation of optical satellite images with ground GPR depth slices or between different geophysical datasets, as well as the impact of noise working with other distributions functions beyond the Gaussian (normal) distribution studied here. This fusion could allow researchers to explore new ways of analysis of existing GPR results based on the continuing increase of the spatial resolution of satellite images and ground-based prospection techniques. A major step towards this direction will be the investigation of potential correlation of satellite imagery (including archive datasets) with large scale geophysical databases. In this step, the impact of scale (e.g., pixel resolution of the satellite image) needs to also be addressed.  [47] p NIR is the near infrared reflectance; p red is the red reflectance; p green is the green reflectance; p blue is the blue reflectance; p x is the reflectance at a specific wavelength.  SG (Sum Green Index) mean of reflectance across the 500 nm to 600 nm [38] p NIR is the near infrared reflectance; p red is the red reflectance; p green is the green reflectance; p blue is the blue reflectance; p x is the reflectance at a specific wavelength.