Geophysical data fusion of ground-penetrating radar and magnetic datasets using 2D wavelet transform and singular value decomposition

This work addresses the problem of the lack of perceptibility that geophysical data may have. Data fusion allows us to combine datasets, providing an improved and more informative source of information about structures buried in the ground. After testing different approaches, a strategy was developed using ground-penetrating radar and magnetic datasets collected over the same area. Data collected at the Roman Villa of Pisões (Beja, Portugal), which is a place of easy application of geophysical methods, were used to test the method, but with problems caused by the properties of the soil. The approach was based on processing operations that allow the fusion of images obtained by different equipment widely used in medical imaging for tumor detection and image processing. The goal is to create an improved image with data fusion that has higher quality than the input images, allowing a better understanding of the object of the study. The approach is composed of two stages: pre-processing and data fusion. Pre-processing is applied to enhance the input data. It consists of removing background noise through singular value decomposition applied in the spectral domain. Then the calculation of the data entropy will highlight the differences corresponding to the spatial alignments compatible with buried structures. Then, both entropy maps of the two datasets are fused in the second processing step to produce the final image. This step involves applying the 2D wavelet transform to each entropy map, decomposing them into sub-bands. Algorithms to calculate multiresolution singular value decomposition and the image gradient are applied to the sub-bands. The processed sub-band pairs are then fused using specific fusion rules. The fused image is obtained by applying the inverse of the wavelet transform. Data fusion with the proposed approach allows us to obtain a detailed image that is sharper and of better quality than the input datasets. The increase in sharpness and quality can be quantified through the sharpness index and the BRISQUE quality index in several steps of the processing. The obtained values confirm the graphical results. Images produced by the proposed data fusion approach suggest that the perceptibility has increased, allowing us to provide conclusions about the existence of buried structures.

1 Introduction

Geophysical surveys implemented in an archaeological environment intend to inspect the subsurface to assess the existence of buried structures. In quite situations, we realize that there are some conditions on the site that prevent suitable detection of these structures. In the ground-penetrating radar (GPR) method, heterogeneities of the soil are given by the soil granulometry, collapses of ancient structures, stones, and obstacles at the surface, such as bushes and holes, that add extra information to the data beyond that corresponding to buried structures such as walls or floors. This additional information is considered noise as it can make it difficult to interpret the results. In addition, in the acquisition, the method subsamples in the direction perpendicular to the profiles due to the distance between the profiles, which is large when compared to the trace separation in a B-scan. In the same way, the magnetic (MAG) method is disturbed by pottery fragments that many times are in the surface. Its existence is an indicator of buried structures. However, they mask the signal due to amplitude peaks.

In several works, despite the strong evidence at the surface, in fact, the results from GPR and MAG surveys do not allow us to clearly understand the existence of buried structures. The lack of perceptibility due to the physical and chemical conditions of the soil and archaeological structures decreases the signal-to-noise ratio of the data. The lack of contrast in the results prevents an effective interpretation of the soil content. In these situations, despite the large amount of noise, we have evidence that the signal has useful information about the buried structures (Oliveira, 2020).

Data fusion is a concept that allows one to perform a joint analysis of two datasets by creating a new dataset that combines the input datasets. The methodology makes it possible to increase the perceptibility of the existence of buried structures. This approach is widely used in medical imaging to detect tumors using different datasets as input.

This work presents an approach of data fusion applied to geophysical data, using GPR and MAG datasets from the same site (spatially overlapped) as input. The data were collected in the Roman Villa of Pisões (Beja, Portugal). The proposed fusion approach uses the 2D wavelet transform, multi-resolution singular value decomposition, and image gradient. The fused data show an increase in sharpness and quality when compared to the input data individually.

2 Roman Villa of Pisões (Beja, Portugal)

2.1 Framework

The Roman Villa of Pisões (Figure 1), located near Penedo Gordo (10 km west Beja), was discovered in 1970 during agricultural works. It is currently owned by the University of Évora and integrated into an experimental farm. This site has the right conditions to perform any geophysical surveys, such as open spaces, without obstacles, and several examples of archaeological structures exposed than can be studied and analyzed in situ (Caldeira and Borges, 2017). Since 2017, several geophysical surveys were performed with different methods to find buried structures in the unburied sectors and to evaluate the state of conservation of the exposed ones.

FIGURE 1

FIGURE 1. Geographical location and general view of the Roman Villa of Pisões (image obtained by an unmanned aerial vehicle). The exposed structures are the remains of the pars urbana of the site.

The structures discovered since 1970 have great opulence (Figure 2) similar to that we can observe in structures such as the big pool, thermal building, mosaic and marble pavements, and walls, defining a large housing complex (Costa, 1983; Alarcão, 1988; Wrench, 1998; Reis, 2004; Couto, 2007; Pereira et al., 2013).

FIGURE 2

FIGURE 2. Several elements that compose the pars urbana of Pisões. (A) Pool with large dimensions. (B) Mosaic pavement with a complex design. (C) Thermal building. (D) Several rooms with walls and pavements.

The splendor and great extent of the unburied structures implies that there were support structures such as the pars rustica (residence for workers) and the pars fructuaria (blacksmithing, pottery, and grainer). These structures are still to be discovered, and knowledge about the unearthed remains is very scarce.

2.2 Geophysical prospection carried out in Pisões

The geophysical study under development in Pisões began in 2017. The first experimental surveys were carried out to gain knowledge about the unknown part of the archaeological site and to create a planification of the actions to protect and conserve the site. The magnetic survey (vertical gradient mode) was applied in several sectors around the pars urbana (Figure 3), to discover the sites of greatest archaeological interest. The GPR survey was applied to only one of the magnetic sectors.

FIGURE 3

FIGURE 3. Location of the prospected sectors. The GPR survey was carried out in sector A1, and the magnetic survey in all sectors.

2.1.1 Magnetic survey

The magnetic method was used to study most of the unexcavated part, in a vertical gradient mode. It was performed in parallel profiles with a spacing of 0.5 m, with control of marks every 10 m.

The team used an array of proton precession overhauser magnetometers to acquire magnetic data in the vertical pseudo-gradient mode (GEMSYS GSM-19). This survey was conducted by scanning all sectors in parallel (zigzag) profiles spaced 0.5 m apart, with markers every 10 m for later lag correction.

Before processing, the raw data had to be analyzed, the operation flow had to be planned, and the parameters for the algorithms had to be determined. However, noteworthy features in the first approximation to the magnetic data are illustrated in Figure 4A: 1) the average values of the magnetic gradient of the anomalies are dispersed over multiple sections of the area, with two ranges of values: one with low values, between ±10 nT/m, and another with high average values. 2) Dependence on the direction of acquiring the profiles that produce a striped pattern on the maps with the same direction of the profile acquisition.

FIGURE 4

FIGURE 4. Processing results of several steps in area A1. The raw data (A) allow one to define the range of amplitude values to be considered. We can observe that most amplitudes are defined between ±20 nT/m (B). The filters mitigated the striped effect (C), and the reduction to the magnetic pole corrects the position of the anomaly map (D).

The processing consists of the application of operations such as the correction of parallax effects, calculation of regular grids, clipping the high values, moving average filter, Gaussian filter, band-pass filter, and reduction to the pole (detailed in Table 1).

TABLE 1

TABLE 1. Processing chain applied to the magnetic data.

First, a study was carried out to define the range of the most significant values for each anomaly map (Figure 4A). Each range of values corresponds to a type of magnetism (Smekalova et al., 2008; Teixidó & Peña, 2018; Oliveira, 2020). The objective is to eliminate all the strongest anomalies attributed to ferric magnetism (> ± 50–100 nT/m) to highlight the induced magnetism anomalies. The most significant values are defined between ±20 nT/m (Figure 4B), corresponding to remanent magnetism (> ± 20–50 nT/m), such as baked clay materials, and induced magnetism, such as masonry structures (< ± 10 nT/m). This step was performed with the software Golden Surfer.

After that, we had to eliminate the background striped pattern with a Gaussian filter and a band-pass filter applied in the spectral domain through the software Geosoft Oasis montaj (Figure 4C). In addition, the reduction to the magnetic pole operation was also applied to the resulting map to verticalize the magnetic anomalies (Figure 4D).

The final anomaly map of area A1 is represented in Figure 4D (Filtered + RTP), where spatial variations in amplitude can be observed that define some alignments that seem to suggest that they are remains of ducts. The magnetic method does not allow matching the data to a depth of investigation. The magnetic data correspond to an integration thickness, whose values are caused by the structures that exist close to the surface.

Applying this processing to all sectors, we obtain the anomaly map of the entire prospected area with the magnetic method. Overlaying the results to an aerial orthophoto, we can interpret all the alignments in the archaeological space, considering the excavated remains of the pars urbana sector (Figure 5).

FIGURE 5

FIGURE 5. Results of vertical gradient magnetic surveys carried out in the seven prospected study areas. All datasets were filtered with identical sequences of operations, adapted to each dataset, with the additional application of reduction to the magnetic pole.

Even with the best processing applied to the magnetic data, we notice a pattern that makes their analysis difficult, which is not very easy to be interpreted if there are buried structures at the site. This is caused by the higher values of magnetic anomaly (spikes), produced by the ceramic fragments that exist scattered throughout the surface and surrounding the structures. Despite this, it is still possible to make an interpretation of the results. The interpretation of the magnetic results will be performed together with the GPR and data fusion results.

2.1.2 Ground-penetrating radar survey

The GPR survey was carried out experimentally only in sector A1, in the same location of the magnetic survey (Figure 3). The choice of this area is due to the magnetic results from A1 sector, which shows two types of anomaly zones: high values (NE part) and low values (SW part). The GPR survey was conducted using a profile spacing of 0.5 m, in a zigzag mode, with the same orientation of the magnetic survey. Due to the apparent lack of quality of the GPR data verified in Pisões, we used two antennas with different central frequency values to evaluate the detection capacity: a 400 MHz antenna (depth range detection between 0.2 and 2.5 m) and a 200 MHz antenna (depth range detection between 0.3 and 4 m). The acquisition parameters for both surveys are detailed in Table 2.

TABLE 2

TABLE 2. Acquisition parameters for GPR surveys.

Before applying the best processing to the GPR data, an analysis of some crucial aspects was performed, such as the signal-to-noise ratio and the frequency range of the reflected signal from the raw data. Figure 6 shows two B-scans of sector A1 obtained in the same location with the two antennas of 400 and 200 MHz, and the corresponding frequency spectrum of each B-scan.

FIGURE 6

FIGURE 6. (A) 400 MHz B-scan GPR-41, without any reflections visible, and the corresponding frequency spectrum, showing several layers of main frequency values, without the central frequency individualized. (B) 200 MHz B-scan GPR-41 (at the same location of the previous one), showing a hyperbolic reflection is very subtle, and the corresponding frequency spectrum, with the central frequency individualized from other values.

The analysis of the two spectra allows us to observe that the reflected energy conserves the radiation pattern of the emitted energy (frequency interval between the values defined by the filters in the acquisition). It was also noted that in the spectrum of the 400 MHz antenna, mono-frequency bands are detected that correspond to multiple reflections caused by the echoes of the air wave (first harmonics of the background at 180, 260, and 420 MHz). Most of the energy is concentrated in these bands, which can explain the lack of useful signal. On the other hand, for the 200 MHz antenna, there is a dense band of reflected energy between 50 and 350 MHz superimposed on the background. It is also observed that the data are centered in the expected band for this antenna, which could indicate the existence of useful signal.

The non-detection of a useful signal observed in the 400 MHz antenna (Figure 6A) may be due to the high conductivity of the Pisões soil (alluvial soils and black clays), which produces a barrier to the propagation of the electromagnetic wave (EMW), in which there is a strong attenuation (Oliveira, 2020). This was confirmed by the calculation of the parameter skin depth (δ). For an attenuation of 1.81 (calculated for this soil considering its wet clay composition: σ = 50 mS/m and ε_r = 27), δ = 0.55 m. At this shallow penetration depth, the estimated percentage of reflected energy is 5.15% of the emitted energy, which means that the medium drastically attenuates the emitted signal and that the energy does not reach the depth of the archaeological stratum, which as will be seen later; its ceiling is located about 0.5 m from the surface.

On the other hand, if we consider an EMW propagation speed of 0.058 m/ns (experimentally determined from the analysis of some B-scans obtained in this area), the central wavelength for the 400 MHz antenna in this soil is λ = 0.14 m. Applying the Huygens criterion, λ/4 = 0.035 m, which corresponds to the size of the particles that can be detected. The average diameter of alluvial soils and silts varies between 0.004 and 0.064 m, so most of these particles will interact with the EMW energy package, producing point reflections in all directions (scattering), causing energy dispersion.

When performing the same calculations for the 200 MHz antenna, a central wavelength of λ = 0.29 m is obtained, with λ/4 = 0.07 m so that less scattering is produced and therefore is more penetrative. This means that with this antenna, it is possible to detect, in this type of terrain, buried structures between 1.0 and 2.5 m deep, as is the case of the GPR-41 profile in Figure 6B. In this B-scan, we can observe a reflection that can correspond to an empty pipe, perhaps a hydraulic structure. The top of the structure is located at about 0.5 m depth.

The GPR data followed four stages of processing:

1) Standard processing (Jol, 2009): 3D-GPR PRC dataset.

2) Standard processing + directional filter (Oliveira, 2020): 3D-GPR1 dataset.

3) Data densification with Fourier interpolation (Oliveira et al., 2022) + standard processing: 3D-GPR2.

4) 3D-GPR2 + noise removal in the spectral domain (Oliveira et al., 2021): 3D-GPR3.

The standard processing comprises operations such as correct position, aerial wave removal, deconvolution, temporal FIR filter, spatial FIR filter, gain adjustment, and Hilbert transform. From the processed 3D-GPR dataset, several depth slices are extracted that can be used as input in the calculation of cover surface (Peña & Teixidó, 2013; Oliveira, 2020).

The cover surface allows the representation of the most expressive reflections compared with the background value of each slice, highlighting the relevant information of all the considered depth slices. This provides a three-dimensional perspective of the reflection alignments. As the cover surface consists of a single slice that integrates information from several depths, it can be a dataset equivalent to the one produced by the magnetic method; that is, both concern an integration thickness that collects information from several levels of depth.

The results of all the stages of the processing are represented in Figure 7. Interpretative lines of the most significant alignments of reflections have been added to aid the interpretation. The interpretive lines are drawn manually through the visual analysis of the results. If there are reflection/anomaly alignments that can define the spatial distribution of a structure, then it could be considered an important reflection. The analysis is carried out by a team of geophysicists who work in an archaeological environment and who have already had evidence of the observations after excavation, in association with a team of archaeologists. In this way, it is possible to establish a parallelism between identical situations.

FIGURE 7

FIGURE 7. Cover surfaces of the results obtained by the advanced processing of the 3D-GPR datasets of sector A1 of Pisões. (A) 3D-GPR PRC, resulting from standard processing. (B) 3D-GPR1, resulting from application of the directional filter. (C) 3D-GPR2, resulting from data densification. (D) 3D-GPR3, resulting from noise removal in the spectral domain. Interpretive lines were drawn on all surfaces, from the most significant alignments, to facilitate their interpretation.

3 Data fusion

3.1 Data fusion concepts

Image fusion has been a subject that has been extensively explored in recent years in several fields of application, such as image analysis and improvement, computer vision (Gautam & Kumar, 2015), and, more relevantly, medical imaging. In this area, data fusion will improve images from complementary diagnostic exams that, by themselves, may not be sufficiently detailed to diagnose the pathology with the necessary rigor and accuracy. The example of image fusion in medical imaging facilitates the understanding of its primary goal. The integration of all the important information in two representations of the same object of study, made from data captured by different examination methods, creates a fused image of the same object that is more informative than the two that gave rise to it and, consequently, with a superior graphic appearance. The third image allows us to achieve the goal of the approach, that is, to provide a better interpretation of the object of study. In medical imaging, this technique can make a difference in the identification and diagnosis of a tumor since it may only be possible to identify it in the image obtained by fusion, making it impossible to be detected in the initial images obtained by different methods.

In the example related in Gautam and Kumar (2015), the images used as input data in the fusion algorithms come from different sensors, namely, the computational tomography (CT) method and the magnetic resonance imaging (MRI) method. CT data highlight information about bones, blood vessels, and soft tissues, while MRI data highlight soft tissue in greater detail. Individually, the results of each of the methods may not provide enough information about the part of the body analyzed that contains bones and tissues. In this way, combining CT and MRI results through data fusion algorithms allows the production of a fused image with all the important information coming from the two initial images. This type of data fusion is called multi-sensor image fusion.

Image fusion algorithms can be classified into three categories (Gautam & Kumar, 2015): pixel-level fusion, fusion in terms of aspects/characteristics, and decision-level fusion. Pixel-level fusion is defined by combining relevant information from each of the images used as input data, pixel by pixel, to generate a composite image that contains more detailed information than the individual input images. This type of fusion is easy to implement numerically, consuming few computational and temporal resources, and is therefore the most used type (Mitianoudis & Stathaki, 2007). Aspect/feature-based fusion uses input image attributes such as color, borders, and textures to select the parts of the two images that will contribute to creating the fused image (Sasikala & Kumaravel, 2007). Decision-level fusion is a type of high-level fusion that combines results from several algorithms to create a final decision on the fusion process of the input images (Tao & Veldhuis, 2009).

In the numerical implementation of the fusion algorithms, these can be applied in the spatial and transformed domains. In the spatial domain, pixel-to-pixel fusing is performed on all input images to obtain a single fused image. The following fusion techniques with implementation in the spatial domain stand out: averaging (Pu, 2000; Pajares & de la Cruz, 2004), principal component analysis (Naidu & Raol, 2008; Sadhasivam et al., 2011), intensity-hue-saturation (Koutsias et al., 2000; Wang et al., 2005), and high-pass filter (Chavez et al., 1991). However, spatial methods sometimes negatively affect the fused images, namely, deformation and reduced contrast. We can overcome this type of problem using fusion approaches that were conceived in the transformed domain. These processes involve the decomposition of images, so this decomposition provides directional information throughout the informative part of the image. Thus, the element of the decomposition contains unique information with different resolutions. The most used mathematical transforms in this type of fusion are the discrete wavelet transform (Li & Manjunath, 1994), stationary wavelet transform (Li et al., 2012), curvelet transform (Ali et al., 2010), and non-subsampled contourlet transform (Wang et al., 2004).

The concept of data fusion was previously applied in geophysics using some of the approaches mentioned, highlighting the work developed by Hilbert et al. (2012) and Karamitrou et al. (2019).

The approach used to implement the data fusion of this work, based on the work of Gautam and Kumar (2015), throughout several stages, combines several techniques for image analysis, namely, the 2D wavelet transform and principal component analysis (Naidu & Raol, 2008). First, the input images are separately decomposed, by calculating the 2D wavelet transform, into low- and high-frequency sub-bands; then, the relevant information for each of them is selected through principal component analysis, and finally, the fusion is performed. The numerical implementation that allows analyzing the principal components of the data selection is carried out using the multi-resolution singular value decomposition technique (Naidu, 2011), which is applied to the low-frequency sub-band. To obtain this, we need to apply the 2D wavelet transform to the input data, which is the approximation part with details such as edges, textures, limits, and changes in image sharpness (Naidu, 2011; Gautam & Kumar, 2015).

Next, we will introduce some of the concepts used in the presented geophysical data fusion approach: discrete wavelet transform, multi-resolution singular value decomposition, the gradient of an image, and fusion rules.

3.2 Discrete wavelet transform

Wavelet transform emerged after 1980 as an alternative to the short-term Fourier transform (Mallat, 1989; Naidu & Raol, 2008). Fourier theory decomposes a signal into a sum of sines and cosines, giving good resolution in the frequency domain. In contrast, in the wavelet theory, the signal is represented by a set of wavelet functions, which offers a good solution in both the time and frequency domains. The wavelet transform is widely used in image processing as it allows the multi-resolution decomposition of an image on a biorthogonal basis, resulting in a non-redundant image representation (Naidu & Raol, 2008).

The application of the discrete wavelet transform to an image I(x,y) consists of the application, in the horizontal direction, of a low-pass filter L and a high-pass filter H, creating the coefficient matrices I _L (x,y) and I _H (x,y) (Gautam and Kumar, 2015) (Figure 11). In the vertical direction, a low-pass filter and a high-pass filter create the sub-bands (sub-images) I _LL (x,y), I _LH (x,y), I _HL (x,y), and I _HH (x,y) (Figure 8). The I _LL (x,y) sub-band represents the approximation part of the image, which contains the average information of the image corresponding to low frequencies. It can be considered a smoothed and subsampled image of the initial image. The sub-bands I _LH (x,y), I _HL (x,y), and I _HH (x,y) represent the detailed part of the input image, containing the directional information in the horizontal, vertical, and diagonal directions, corresponding to the high frequencies. Multi-resolution is obtained if the wavelet transform is successively applied to the low-frequency sub-band. The decomposed image can be reconstructed by applying the inverse of the wavelet transform.

FIGURE 8

FIGURE 8. Representative scheme of the application of the 2D wavelet transform to an image (Gautam and Kumar, 2015), in which the following are shown: (A) one level of decomposition and (B) two levels of decomposition.

3.3 Multi-resolution singular value decomposition

Multi-resolution singular value decomposition (MSVD) is a case of implementing the SVD technique to data in which the 2D wavelet transform was applied. In wavelet transforms, the signal from the decomposition sub-bands is filtered separately with low-pass and high-pass filters. In the case of MSVD application, instead of applying filters to the sub-bands, the SVD technique is used (Naidu, 2011) to select the useful information of each sub-band.

When applying MSVD to a sub-band, two matrices are obtained: a matrix U, which contains the singular vectors of the sub-band, and another matrix, where each row defines, respectively, and in order, the approximation part Φ of the sub-band (first row), corresponding to the highest singular values, and the detailed parts Ψ (next three rows), corresponding to the remaining singular values (Figure 9). The decomposed image can be reconstructed by applying the inverse of the MSVD.

FIGURE 9

FIGURE 9. MSVD decomposition structure of an image with three levels of decomposition (Naidu, 2011).

3.4 Image gradient

The gradient calculation was designed to fuse images whose goal is the selection of the sharpest pixels from the sets of images to be fused. This allows us to build a new image with higher quality, with pixels better focused (Yang et al., 2014). The authors of this approach conceived a measure of focus of sharpness that allows selecting the desired coefficients in the matrices that define the sub-bands obtained by applying the discrete wavelet transform. The sharpness measure is a Tenengrad function based on a Sobel operator that combines information from neighboring pixels through a fixed window that traverses the entire data matrix (Yang et al., 2014). The Sobel operator on which the Tenengrad function is based refers to the Prewitt operator, which is used in signal processing to detect edges and contours of an image (Chaple et al., 2015). Technically, it is a discrete differentiation operator that calculates an approximation of the gradient of the image intensity function in the horizontal and vertical directions. This allows one to determine the direction of the greatest possible increase from light to dark color and the respective rate of change. The sharpening measure can select the pixels in an image that give it the most sharpness.

The fusion method, through gradient calculation proposed by Yang et al. (2014), is used to fuse the sub-bands corresponding to the high frequencies obtained by applying the discrete wavelet transform to two images (LH, HL, and HH sub-bands). The numerical implementation requires the calculation of the horizontal and vertical gradient coefficients for each sub-band of the two images. Each gradient coefficient corresponds to the directional derivative in the horizontal and vertical directions, respectively. The MATLAB function getGradientH, from Paul, Sevcenco, and Agathoklis (2016), allows one to calculate it. Then, the calculated coefficients are combined by calculating the magnitude of the gradient vector.

3.5 Fusion rules

Fusion rules consist of applying mathematical operations whose result will be the matrix that has the fused data. These consist of simple operations, such as calculating the average or maximum matrix, pixel by pixel. Depending on the rule used, the result will be different. For example, if we consider the matrix with the sharpest data from each dataset, the selection of the sharpest data from both datasets can be done by calculating the maximum of each pixel of both matrices.

3.6 Data fusion of ground-penetrating radar e magnetic datasets

The proposed approach to fuse GPR and MAG datasets was based on the work developed by Gautam and Kumar (2015), which describes a methodology for image fusion obtained by different equipment, combining 1) 2D wavelet transform, 2) multi-resolution singular value decomposition, and 3) the gradient calculation.

The proposed numerical scheme consists of two stages: pre-processing and data fusion. The need to incorporate a pre-processing step is because the input datasets, even in the situation of high perceptibility, are not capable of providing results that can show the effectiveness of the proposed approach.

The following datasets were used as inputs: cover surface obtained by the GPR method using the 200 MHz antenna and a map of magnetic anomalies obtained by the MAG method. Both datasets correspond to an integration thickness of several levels deep and are therefore considered equivalent to be used as inputs to the data fusion approach.

3.7 Pre-processing

The sequence of operations that constitute the pre-processing operation, applied to each dataset, is composed of the following procedures:

1) Application of a median filter to smooth the datasets.

2) Automatic alignment of the two datasets to ensure they are precisely overlapping.

3) Applying the SVD technique in the 2D Fourier domain to select only the useful information, excluding background noise (Oliveira et al., 2021).

4) Calculation of the local entropy of the data so that only the amplitudes corresponding to buried structures are highlighted. The 2D data entropy shows the spatial disorder of the data. In the GPR and MAG data, the background corresponds to the organized part and the buried structures to the disordered part.

3.8 Data fusion

The sequence of operations that allows the fusion of geophysical datasets is summarized below and is outlined in Figure 10:

1) Application of the 2D wavelet transform, obtaining the following sub-bands: LL, HL, LH, and HH.

2) The MSVD algorithm is applied to the LL sub-band, obtaining the following matrices: U, Φ, Ψ₁, Ψ₂, and Ψ₃.

3) The fusion operations are carried out as follows:

1) The matrices U and Φ, from each dataset, are fused by calculating the average matrix, pixel by pixel, obtaining fused matrices U_F and Φ_F with the exact dimensions as the unfused matrices.

2) The matrices Ψ₁, Ψ₂, and Ψ₃, from each dataset, are fused by calculating the maximum matrix, pixel by pixel, resulting in the fused matrices Ψ_1F, Ψ_2F, and Ψ_3F, with the exact dimensions as the unfused matrices.

3) The fused matrices U_F, Φ_F, Ψ_1F, Ψ_2F, and Ψ_3F are used by the algorithm that applies the inverse of the MSVD function to calculate the LLF matrix, with the exact dimensions as the unfused matrices.

4) The LH, HL, and HH sub-bands of each dataset are fused by calculating the gradient of each sub-band of each dataset. The fusing operation is performed by calculating the average matrix, pixel by pixel, obtaining the fused LH_F, HL_F, and HH_F matrices, with the exact dimensions as the unfused matrices.

5) The fused matrices LL_F, LH_F, HL_F, and HH_F are used by the function that applies the inverse of the wavelet transform, obtaining the final matrix of the fused data with the exact dimensions as the unfused matrices.

FIGURE 10

FIGURE 10. Schematic of the numerical implementation of the geophysical data fusion of GPR and magnetic datasets. It was adapted from Gautam and Kumar (2015) for this fusion approach.

3.9 Evaluation of the results

The evaluation of the results of this methodology can be done by establishing a graphical comparison between the datasets and comparing the values of parameters such as the sharpness index and the BRISQUE quality index at several steps of the processing.

The sharpness index (SI) quantifies the sharpness of an image with a dimension (m × n). To calculate it, we can use the Hudgin gradient (dH) of the image (Paul et al., 2016), which is used in image fusion to select the sharpest pixels (Yang et al., 2014), considering the neighboring pixels. The numerical implementation requires calculating the horizontal components of the Hudgin gradient dxH and vertical dyH. These coefficients are needed to calculate the gradient modulus, used in the SI calculation (Birdal, 2022), through Eq. 1.

$Sharpnessindex=100\ast {\displaystyle \sum }_{i,j}^{m,n}\frac{{\displaystyle \sum }_{i,j}^{m,n}\left|{dH}_{i,j}\right|}{m\times n},(1)$

were m × n is the dimension of the matrix, and i,j define the coordinates of each cell.

The increase in the Hudgin gradient values of the image occurs the more contrasted the image is. We can consider that the sharpness of an image must be associated with the diversity of contrasts. The more significant and more abundant the contrasts, the higher this index.

The BRISQUE quality index allows one to calculate a value, between 0 and 100, for the absolute quality without any reference, using as an evaluator a database with natural scenes with distortions, with artifacts, unsharpness, and noise (Mittal et al., 2011, 2012). The lower value corresponds to the best quality.

4 Results

4.1 Fusion results

The proposed approach to perform geophysical data fusion using different datasets from GPR and MAG methods will be evaluated by the analysis of the inputs and outputs.

The inputs were processed according to standard processing, using the best parametrization to produce the better results. As the input data, we used the 3D-GPR3 dataset (Figure 11A) and the processed MAG dataset (Figure 11B), both acquired in sector A1. Figure 11 shows the input data for the algorithm with normalized amplitude values.

FIGURE 11

FIGURE 11. Input data of the geophysical data fusion algorithm: (A) 3D-GPR3—Coverage surface of the densified, processed, and filtered data from the background noise and (B) filtered vertical gradient MAG with RTP application.

The result obtained with the fusion approach is presented in Figure 12A. Unfortunately, the color scale used in the representation of the normalized gradient colormap presents a slight lack of contrast, which may make it challenging to identify the alignments of possible structures.

FIGURE 12

FIGURE 12. (A) Result of data fusion. (B) Histogram of the absolute frequency of the amplitude values of the results. (C) Adjustment of the color scale of the data fusion result—representation of 95% of the values. (D) Histogram of the absolute frequency of the amplitude values of 95% of the values considered in (A).

The histogram of the amplitude values was calculated to help equalize the color scale to optimize the interpretation of the data without harming the signal itself (Figure 12B). If 95% of the data are considered (mean ±2 standard deviation), the final map has a color scale with more contrast, facilitating its interpretation (Figure 12C). The histogram for 95% of the data maintains the initial distribution of the absolute frequency of the amplitude values (Figure 12D). The two histograms present a well-defined normal distribution close to the calculated theoretical model (red line).

The results were evaluated at several processing steps. The sharpness index (Table 3 and Figure 13A) increased, meaning there was an increase in sharpness with the application of the fusion algorithm. On the other hand, the BRISQUE quality index (Table 3; Figure 13B) decreased, meaning that the quality increased with the application of the fusion algorithm.

TABLE 3

TABLE 3. Sharpness index and BRISQUE quality index of sector A1 data at several processing stages.

FIGURE 13

FIGURE 13. Graphical representation of the evolution of the sharpness index (A) and the BRISQUE quality index (B) of the data considered in the several stages of the data fusion.

4.2 Interpretation of the results

Figure 14 shows the geophysical models obtained for sector A1: 1) magnetic anomalies, 2) 3D-GPR3 dataset, and 3) fused dataset. The procedure to interpret the geophysical anomalies consists of drawing the most significant alignments of the two input images, as well as the same for the fused dataset. After that, we can make the comparison between the interpretations and make the final interpretation with all possible sources of information.

FIGURE 14

FIGURE 14. Visualization of the models obtained in sector A1. (A) Magnetic anomalies. (B) 3D-GPR3 dataset. (C) Fused data. We drew the most significant interpretive lines on the three models. The fused dataset has more information available than the inputs individually.

In sector A1, we can make the following interpretation, by analyzing Figure 14C:

1) E1: The combination of the hyperbolic alignment visible in the GPR results, the high values of magnetic anomaly and its geometric shape, suggests that it can correspond to a brick wall with a thickness between 1 and 2 m. This structure can be observed in Figure 6B. The top of the structure is located at about 0.5 m depth.

2) E2: This structure appears to be destroyed, so unsure about the interpretation. It has a circular shape, with an inner diameter of approximately 30 m and a thickness of about 5 m (in yellow). These alignments are not compatible with the typical geometry of the Roman structures. The NE half appears to be better preserved; however, some material can be absence that could indicate subsequent reuse. In this case, the high magnetic anomalies could indicate accumulations of bricks and places of combustion.

3) E3: It may be a set of structures that occupy the lower SE quarter of the study area. Given the size and arrangement of the anomalies, they may be the remains of rectangular enclosures of varying dimensions, most of which do not exceed 4 m × 4 m. In some of them, point magnetic anomalies were detected, which are compatible with housing areas and some remains of a domestic oven since the magnetism is not very high.

4) E4: Inside the circular crown (E2), scattered traces of enclosures of different sizes and typologies were detected. Toward the NE part, the larger and better preserved remains are located, which may correspond to the corners of houses or ovens, while in the rest of the circle, the structures seem to be more destroyed.

The interpretation of sector A1 was carried out from a geophysical point of view that tracks the possible structures that generate anomalies. However, in a complex and polyphasic place such as Pisões, an interpretation using archaeology that considers the urban canon and its chronology is required. Despite this lack, the alignments interpreted in A1 were used as a standard to extrapolate the interpretation of all the magnetic areas.

In Figure 15 a general interpretation is represented in which the models were georeferenced and superimposed on an orthophoto of the archaeological site. This interpretation is approximate and should be considered in its generic context, so many of the structures described will have to be corrected after making the comparison with future excavations.

FIGURE 15

FIGURE 15. Interpretation of the geophysical prospection on the orthophoto of Villa of Pisões.

Considering the values and the spatial distribution of the magnetic anomalies, the archaeological site presents a zone of high anomalies (± 60 nT/m), located in the northern part, and another zone of medium anomalies, located in the south (± 20 nT/m). Separating the two environments, a tripolar wide linear anomaly (blue ditch) was detected, characterized by a central maximum and two lateral minima, which was interpreted as a hydraulic channel in which the maximum corresponds to the cemented part (opus signinum?), with a thickness between 0.5 and 1 m, and the lateral parts to the excavated area between 1 and 1.5 m transversal to each side, later filled.

In sector A2, an indeterminate form appears to begin, which coincides with a slight gradient at altitude 185 m (blue line in Figure 15). On its right margin, a scattered set of high anomalies with an average dimension of 3 m × 4 m was detected, which may correspond to enclosures that contain some type of cementation (opus signinum? or opus caementicium?), while on the left margin, smaller anomalies are located (approximately 2 m × 2 m), possibly formed by bricks, which in addition to the current path seem to correspond to wells. Almost parallel to the channeling, a negative rectangular anomaly was detected that can be associated with a ditch. The south corner is also characterized by a negative anomaly with a similar interpretation, possibly related to the previous wells located at a higher altitude.

In sector A3, a series of anomalies were obtained, most of which must have an edaphological-biological origin.

In sector A4, there is a linear tripolar anomaly that seems to end at the end of sector A5, a kind of enclosure about 5 m wide. As in sector A2, here the channeling appears to be a divide between structures with high magnetic values ​with similar dimensions and constructions, and in the south zone with medium anomalies that extend to sectors A6 and A7. In the structures to the north, small circular bipolar anomalies were located that may correspond to ovens (marked with red circles).

If the structures located on the eastern edge of the villa are considered, with the average anomalies detected in the A6 and A7 sectors, it seems likely that they are similar housing spaces. In sector A7, a similar anomaly was detected that seems to come out of the riverbed toward the great natatio.

5 Discussion

5.1 Geophysical data fusion

In this work, a problem of lack of perceptibility observed was presented in the results of geophysical surveys. The apparent lack of information is related with the low contrast between signal and noise, due to the physical and chemical conditions of the ground. The initial hypothesis of the proposed approach consists of the assumption that the data hold useful information mixed with noise. A possible solution is the combination of datasets from different geophysical methods to create a new dataset with more information than the inputs. Data fusion is a type of approach widely used in medical imaging and image processing to improve the input datasets by creating a new one with better quality.

In the standard processing of GPR and MAG datasets, data fusion is not a common operation. The bibliography of the GPR method shows that it is necessary to apply advanced processing to extract more information from the data. This type of operation uses mathematical tools that allow one to convert the data to transformed domains, facilitating the application of filters due to the symmetry characteristics of the transformed data, such as the case of the 2D Fourier transform. Decomposition techniques, such as singular value decomposition and the 2D wavelet transform, are widely used to decrease the noise and increase the sharpness of images, respectively. This type of operation is applicable to geophysical datasets. We highlight the use of SVD to erase noise in GPR data, applied in the 2D Fourier domain (Oliveira et al., 2021), and now the use of SVD at a multi-resolution level and the wavelet transform to perform data fusion of geophysical datasets.

In the presented case of study of Pisões, the fused dataset shows anomaly alignments that suggest there could be correspondence with buried structures. The fused dataset has more information than the inputs individually. Both GPR and MAG data do not allow conclusions to be drawn. Alignments detected in the fused dataset have preferred directions compatible with direction three identified in previously excavated structures (Figure 16). The structures in direction three do not correspond to the directions of the structures that make up the pars urbana, nor are they the same as those observed in the baths. This may indicate the location of the support structures of the villa that have not yet been discovered but that are assumed to have existed given the splendor of the pars urbana.

FIGURE 16

FIGURE 16. Interpretation of the anomaly alignments of the fused dataset regarding the direction of the excavated structures. The fused alignments had correspondence with the direction 3, different pars urbana directions, and the baths. This could indicate the location of the support structures of the main part of the Roman Villa that have not yet been discovered.

The geophysical fusion approach consists of the application to each used dataset (GPR and MAG) of the 2D wavelet transform that decomposes the input data into sub-bands. To each one is applied MSVD and calculated the gradient. The resulting sub-band pairs are fused with specific fusion rules. The final operation is the application of the inverse of wavelet transform to produce the fused dataset. This approach allows one to obtain a more detailed image, with more information than the inputs. The sharpness and quality indexes calculated prove the results in several steps of the processing.

5.2 Potential limitations

The proposed approach has some limitations, namely, the impossibility of assigning weights to each dataset. For example, in the Pisões results, MAG data have more useful information than GPR. The amount of information from each method is not quantified in this proposed approach. This can be an important step, as the most informative data can carry greater weight than the apparently less informative data. However, this consideration goes against the premise that we considered in the fusion proposal, which says that the data apparently without useful information can contain relevant information that can contribute positively to the outcome of the data fusion. The proposed fusion method considers the weight of each one when applying the maximum function to the SVD sub-band that allows one to select more informative pixels. We think that if there was a way to quantify the useful information in the data to calculate the respective weight of importance, it can increase the effectiveness of the methodology.

Another limitation is the high influence of the input data to produce a fused image with reasonable quality. Currently, it is necessary to amplify the useful information of the input datasets through a pre-processing step that includes the use of SVD and the calculation of the entropy of the data. By applying these operations, although its effectiveness, this can modify the data with the risk of the information becoming corrupted.

Another limitation is the characteristic slowness of an archaeological excavation, which prevents it from being able to quickly prove the results obtained by geophysical surveys. For example, at the Pisões site, the fused results have some weird alignments that need to be verified by an excavation. There is currently no plan for this to happen.

5.3 Contribution of the current understanding

The proposed methodology to apply data fusion to geophysical datasets, from GPR and MAG methods, is a new approach to enhance the data from these methods. Data fusion is widely applied in medical imaging and image processing, and there are several ways to implement it. Application to geophysics is quite new and is not a type of standard operation from the processing software of each method.

Usually, to conclude about structures and discontinuities in the ground, it is necessary to apply several methods to cross all the information provided by the results. The global interpretation is made by the individual analysis of the results of each method. However, sometimes the available information is not sufficient to be able to draw a conclusion about the existence of buried archaeological structures. Data fusion is an interesting concept that can help the global interpretation of the results of several methods applied in the same location. The proposed approach allowed implementing data fusion to GPR and MAG datasets, and a new image obtained from both inputs was created, with more information than the inputs, and with more quality and sharpness. This allowed us to increase our knowledge about the Pisões site, where both GPR and MAG methods had failed due to the terrain conditions.

5.4 Future direction of the research

In the future work, we intend to extend the fusion methodology to other geophysical methods, such as seismic, electrical resistivity tomography, and electromagnetic induction. We can use other datasets as input pairs like the presented case of study or use several methods at the same time to create a much better result.

We also intend to solve the presented limitations regarding the influence of input data and the weights of each input dataset. As for the slowness of the excavations, we can only contribute with images as realistic as possible so that this can motivate the occurrence of an excavation in a short time.

6 Conclusion

This study provides a new methodology to perform data fusion of geophysical datasets, using ground-penetrating radar (GPR) and magnetic (MAG) data as inputs. Geophysical data fusion allows us to increase the available information regarding structures or discontinuities in the ground. The case of study presented is from an archaeological site, the Roman Villa of Pisões (Beja, Portugal), which is a good place to perform any geophysical survey. GPR and MAG data obtained there apparently did not have useful information to allow them to conclude the existence of buried structures, due to its physical and chemical conditions. However, at the site, there are several clues that suggest the existence of buried remains. We assumed that the information exists mixed with the noise, so the information can be highlighted with processing. Data fusion, which is an approach widely used in medical imaging and image processing, allows the creation of a new dataset using other datasets as input, with more information and quality than the inputs individually. The application to geophysical datasets was effective and allowed us to obtain information about the subsurface of Pisões that was absent before the fusion results.

The geophysical data fusion is implemented using decomposition algorithms, such as 2D wavelet transform and multi-resolution singular value decomposition, as well as image gradient. All the algorithms are applied successively to obtain decomposition of the signal into sub-bands iteratively. The obtained result is a new image, which is sharper and has superior quality than the input datasets. To evaluate the results, the sharpness index and BRISQUE quality index were calculated. These values indicate the increase in sharpness and quality, contributing to the validation of this approach.

The use of data fusion techniques applied to geophysical datasets, using mathematical transforms and automatic data selection, allows us to enhance geophysical datasets and consequently improve the obtained results and interpretations that we can perform. However, we should highlight that, currently, there are some limitations, such as the impossibility to give weights to each input dataset, related to the level of available information in the data, the high influence of the inputs in the production of the fused results, and the slowness of the archaeological excavations that prevent us from proving the geophysical results in a short time. The first two are possible to solve following the investigation in this topic.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

RO led this research with BC and TT. RO, BC, JB, MB and TT contributed to data gathering and modeling. RO contributed to the calculations and computational results. RO, BC, JB, MB and TT contributed to data analysis, discussion, and manuscript preparation and revision.

Funding

This work has been partially supported by the research project “Innovación abierta e inteligente en la EUROACE” with the reference 0049_INNOACE_4_E, by the European Union through the European Regional Development Fund included in COMPETE 2020, and through the Portuguese Foundation for Science and Technology (FCT) by the projects UIDB/04683/2020-ICT (Institute of Earth Sciences) and SFRH/BSAB/143063/2018.

Acknowledgments

The authors express their gratitude to the archaeologists André Carneiro (University of Évora, Portugal) and Jesús García Sánchez (Institute of Archaeology of Mérida, Spain) for their contribution to the interpretation of geophysical anomalies in an archaeological context.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alarcão, J. (1988). Roman Portugal. Introduction I, gazetteer II/fasc. 2. Warminster, England: Coimbra e Lisboa. Aris and Phillips.

Google Scholar

Ali, F. E., El-Dokany, I. M., Saad, A. A., and Abd El-Samie, F. E. (2010). A curvelet transform approach for the fusion of MR and CT images. J. Mod. Opt. 57 (4), 273–286. doi:10.1080/09500340903541056

CrossRef Full Text | Google Scholar

Birdal, T. (2022). MATLAB central file exchange. Available at: https://www.mathworks.com/matlabcentral/fileexchange/32397-sharpness-estimation-from-image-gradients (accessed on April 13, 2022).Sharpness estimation from image gradients

Google Scholar

Caldeira, B., and Borges, J. F. (2017). Relatório de Geofísica Prévio à Construção do Centro de Interpretação e Acolhimento do Sítio Arqueológico da Villa Romana de Pisões. Évora, Portugal: Universidade de Évora.

Google Scholar

Chaple, G. N., Daruwala, R. D., and Gofane, M. S. (2015). Comparisions of Robert, Prewitt, Sobel operator based edge detection methods for real time uses on FPGA. Proc. - Int. Conf. Technol. Sustain. Dev. ICTSD 2015, 1–4. doi:10.1109/ICTSD.2015.7095920

CrossRef Full Text | Google Scholar

Chavez, P. S., Sides, S. C., and Anderson, J. A. (1991). Comparison of three different methods to merge multiresolution and multispectral data. Landsat TM and SPOT panchromatic. Photogrammetric Eng. Remote Sens. 57 (3), 295–303.

Google Scholar

Costa, M. L. V. (1983). Contribuições para o estudo de alguns dos mosaicos da Villa Romana de Pisões. Arq. Beja, 2a Série (2), 95–122.

Google Scholar

Couto, M. B. S. M. (2007). Balnevm da villa romana de Pisões. Lisboa, Portugal: Universidade Nova de Lisboa.

Google Scholar

Gautam, S., and Kumar, M. (2015). An effective image fusion technique based on multiresolution singular value decomposition. INFOCOMP 14 (2), 31–43.

Google Scholar

Hilbert, C., Grandjean, G., Bitri, A., Travellentti, J., and Malet, J. (2012). Characterizing landslides through geophysical data fusion: Example of the La Valette landslide (France). Eng. Geol. 128, 23–29. doi:10.1016/j.enggeo.2011.05.001

CrossRef Full Text | Google Scholar

Jol, H. M. (2009). Ground penetrating radar. Elsevier, 524.

Google Scholar

Karamitrou, A., Bogiatzis, P., and Tsokas, G. (2019). Fusiono f geophysical images in the study of archaeological sites. Archaeol. Prospect. 27, 119–133. doi:10.1002/arp.1766

CrossRef Full Text | Google Scholar

Koutsias, N., Karteris, M., and Chuvieco, E. (2000). The use of Intensity-Hue-Saturation transformation of Landsat-5 Thematic Mapper Data for burned land mapping. Photogrammetric Eng. Remote Sens. 66 (7), 829–839.

Google Scholar

Li, H., and Manjunath, B. s. (1994). Multi-sensor image fusion using the wavelet transform, 11.

Google Scholar

Li, Z., Chai, Y., Guo, M., and Li, H. (2012). Multifocus image fusion scheme based on feature contrast in lifting stationary wavelet domain. Chongqing Daxue Xuebao/Journal Chongqing Univ. 35 (10), 109–116.

Google Scholar

Mallat, S. G. (1989). A theory for multiresolution signal decomposition: The wavelet representation. IEEE Trans. Pattern Anal. Mach. Intell. 11 (7), 674–693. doi:10.1109/34.192463

CrossRef Full Text | Google Scholar

Mitianoudis, N., and Stathaki, T. (2007). Pixel-based and region-based image fusion schemes using ICA bases. Inf. Fusion 8, 131–142. doi:10.1016/j.inffus.2005.09.001

CrossRef Full Text | Google Scholar

Mittal, A., Moorthy, A. K., and Bovik, A. C. (2011). Blind/referenceless image spatial quality evaluator. 2011 conference record of the forty fifth asilomar conference on signals. Syst. Comput. (ASILOMAR) 727, 723. doi:10.1109/ACSSC.2011.6190099

CrossRef Full Text | Google Scholar

Mittal, A., Moorthy, A. K., and Bovik, A. C. (2012). No-reference image quality assessment in the spatial domain. IEEE Trans. Image Process. 21 (12), 4695–4708. doi:10.1109/tip.2012.2214050

PubMed Abstract | CrossRef Full Text | Google Scholar

Naidu, V. P. S. (2011). Image fusion technique using multi-resolution singular value decomposition. Def. Sci. J. 61 (5), 479–484. doi:10.14429/dsj.61.705

CrossRef Full Text | Google Scholar

Naidu, V. P. S., and Raol, J. R. (2008). Pixel-level image fusion using wavelets and principal component analysis. Def. Sci. J. 58 (3), 338–352. doi:10.14429/dsj.58.1653

CrossRef Full Text | Google Scholar

Oliveira, R. J., Caldeira, B., Teixidó, T., Borges, J. F., and Carneiro, A. (2022). Increasing the Lateral Resolution of 3D-GPR Datasets through 2D-FFT Interpolation with Application to a Case Study of the Roman Villa of Horta da Torre (Fronteira, Portugal). Remote Sens. 14 (16), 4069. doi:10.3390/rs14164069

CrossRef Full Text | Google Scholar

Oliveira, R. J., Caldeira, B., Teixidó, T., and Borges, J. F. (2021). GPR clutter reflection noise-filtering through singular value decomposition in the bidimensional spectral domain. Remote Sens. 13 (10), 2005. doi:10.3390/rs13102005

CrossRef Full Text | Google Scholar

Oliveira, R. J. (2020). Prospeção Geofísica aplicada à Arqueologia. PhD Thesis in Earth and Space Sciences (Geophysics). Évora, Portugal: Institute of Research and Advanced Training, University of Évora.

Google Scholar

Pajares, G., and de la Cruz, J. M. (2004). A wavelet-based image fusion tutorial. Pattern Recognit. 37 (9), 1855–1872. doi:10.1016/s0031-3203(04)00103-7

CrossRef Full Text | Google Scholar

Paul, S., Sevcenco, I. S., and Agathoklis, P. (2016). Multi-exposure and multi-focus image fusion in gradient domain. J. Circuits, Syst. Comput. 25 (10), 1650123–3. doi:10.1142/s0218126616501231

CrossRef Full Text | Google Scholar

Peña, J. A., and Teixidó, T. (2013). Archaeological applications. RNM104 - informes, 10. Granada, Spain: Universidad de Granada - Instituto Andaluz de Geofísica, 1.Cover surfaces as a new technique for 3D GPR image enhancement

Google Scholar

Pereira, C., Soares, A., and Soares, R. (2013). Os mausoléus da villa romana de Pisões: A morte no mundo rural romano. Rev. Port. Arqueol. 16, 303–321.

Google Scholar

Pu, T. (2000). Contrast-based image fusion using the discrete wavelet transform. Opt. Eng. 39 (8), 2075. doi:10.1117/1.1303728

CrossRef Full Text | Google Scholar

Reis, M. P. (2004). Las termas y balnea romanos de Lusitania. Madrid: MNAR.

Google Scholar

Sadhasivam, S. K., Keerthivasan, M. B., and Muttan, S. (2011). Implementation of max principle with PCA in image fusion for surveillance and navigation application. Electron. Lett. Comput. Vis. Image Analysis 10 (1), 1–10. doi:10.5565/rev/elcvia.353

CrossRef Full Text | Google Scholar

Sasikala, M., and Kumaravel, N. (2007). A comparative analysis of feature based image fusion methods. Inf. Technol. J. 6, 1224–1230. doi:10.3923/itj.2007.1224.1230

CrossRef Full Text | Google Scholar

Smekalova, T., Voss, O., and Smekalov, S. (2008). Magnetic surveying in Archaeology. More than 10 years of using the Overhauser MSG-19 gradiometer. S. Petersburg: Publishing House of the Polytechnical University.

Google Scholar

Tao, Q., and Veldhuis, R. (2009). Threshold-optimized decision-level fusion and its application to biometrics. Pattern Recognit. 42 (5), 823–836. doi:10.1016/j.patcog.2008.09.036

CrossRef Full Text | Google Scholar

Teixidó, T., and Peña, J. (2018). Ayuntamiento de Málaga. Gerencia municipal de Urbanismo y obras públicas. Ref. AGA-215, informe interno. Área de Geofísica aplicada. España: Instituto Andaluz de Geofísica, Universidad de Granada.Servício para documentación del Yacimiento Fenicio de Cerro del Villar com métodos geofísicps

Google Scholar

Wang, Z., Bovik, A. C., Sheikh, H. R., and Simoncelli, E. P. (2004). Image quality assessment: From error visibility to structural similarity. IEEE Trans. Image Process. 13 (4), 600–612. doi:10.1109/tip.2003.819861

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Zhijun, Ziou, D., Armenakis, C., Li, D., and Li, Q. (2005). A comparative analysis of image fusion methods. IEEE Trans. Geosci. Remote Sens. 43 (6), 1391–1402. doi:10.1109/tgrs.2005.846874

CrossRef Full Text | Google Scholar

Wrench, I. (1998). Transição de signos: Os motivos representados num mosaico da villa romana de Pisões. Rev. Fac. Ciências Sociais Humanas 1 (11), 247–255.

Google Scholar

Yang, Y., Huang, S., Gao, J., and Qian, Z. (2014). Multi-focus image fusion using an effective discrete wavelet transform based algorithm. Meas. Sci. Rev. 14 (2), 102–108. doi:10.2478/msr-2014-0014

CrossRef Full Text | Google Scholar