The Menkaure Pyramid is the smallest of the three main pyramids on the Giza Plateau. Recently, the possibility of a second entrance to the Pyramid has been hypothesized by Van den Hoven [1], based on similarities between the polished granite blocks covering the Eastern face and the blocks around the main entrance on the Northern face. To test this hypothesis, measurement campaigns using three non-destructive techniques, Electrical Resistivity Tomography (ERT), Ground Penetrating Radar (GPR), and Ultrasonic Testing (UST), were carried out on the Eastern face of Menkaure Pyramid. ERT data was obtained from measurements of four long parallel profiles using stainless steel mesh electrodes and inverted using a three-dimensional inversion algorithm. These ERT results guided the more focused grid measurements of a dual-frequency GPR instrument (200/600 MHz antenna) and a 16-channel UST array device. Image Fusion (IF) was utilized to merge the reconstructed ERT, GPR, and UST images, allowing for precise correlation of the detected features from each technique. The images revealed two anomalies directly behind the polished granite blocks, which would indicate the presence of air-filled voids. This interpretation was supported by a series of numerical simulations that considered various possible scenarios under real-world conditions.
Keywords
ScanPyramids
Menkaure pyramid
Non-destructive testing (NDT)
Electrical resistivity tomography (ERT)
Ground penetrating radar (GPR)
Ultrasonic testing (UST)
Image fusion (IF)
1. Introduction
The Menkaure Pyramid is the smallest of the three main pyramids on the Giza Plateau and was named "Menkaure is Divine" [1] (Fig. 1a). Its pyramid complex has been the least disturbed since the Old Kingdom. The builders left the stone pyramids unfinished, and the temples were completed in mud brick. From 1906 to 1910, Reisner excavated the pyramid complex and uncovered evidence revealing a more intricate and dynamic history than its architectural design suggests. This archaeological insight into the life cycle of a pyramid has raised as many questions — if not more — than those associated with the two larger pyramids at Giza. Archaeologists and scientists have always dedicated their investigations to the Khufu Pyramid, but no work has been done on the Menkaure Pyramid since Reisner's work. Therefore, the ScanPyramids team must open a new investigation to reveal more secrets about this pyramid.
Fig. 1. (a) the Northern face of the Menkaure Pyramid and (b) a close-up view of the main entrance, and (c) the polished part of the Eastern face (subject of interest in this study).
One of the defining elements in the original design of the Menkaure Pyramid is its granite cladding. The design initially aimed to cover the entire structure, but it was ultimately limited to only 16 to 18 horizontal courses of stone blocks. Today, just seven rows remain preserved at the base of the pyramid. Recently, the potential existence of a second entrance to the pyramid has been proposed based on visual similarities between sections of granite cladding on the Eastern face and those around the main entrance on the Northern face (Fig. 1b–c) [2]. While most granite blocks have roughly finished surfaces, these two areas stand out with well-polished, accurately fitted blocks. To investigate the potential presence of void spaces behind the granite blocks on the Eastern face of the Menkaure Pyramid, which could support the hypothesis of a second entrance, three non-destructive techniques (NDT) were used: Electrical Resistivity Tomography (ERT), Ground Penetrating Radar (GPR), and Ultrasonic Testing (UST).
A combination of NDT methods is often used in archaeological research to overcome the physical limitations of each method in terms of resolution, depth of investigation, and coverage, and to validate interpretations. Recently, an unknown void hidden behind the Chevron structure on the Northern face of the Khufu Pyramid was detected using various muon techniques [3]. Subsequently, through further detailed and focused muography measurements, the void was determined to be a 9-m-long corridor, identified as the ScanPyramids North Face Corridor (SP-NFC) [4]. The combination of GPR and UST measurements in conjunction with image fusion (IF) provided an estimation of the corridor's shape, cross-sectional dimensions, and precise location [5]. Ultimately, the integrated approach successfully defined the detected void as a corridor with a chevron-shaped ceiling, measuring approximately 2.30 m in height and 2.15 m in width [5].
The ERT method is typically applied in the initial stage of archaeological investigations to quickly identify anomalies that require further detailed study [6]. For example, this approach was used to search for ancient artifacts buried on the western bank of the Nile Valley in Luxor, Egypt [7]. An area of 140 m × 60 m was initially surveyed using ERT and shallow seismic refraction methods. As a result, the survey area was narrowed down to 20 m × 40 m, where GPR was used to accurately determine the location of a statue with dimensions of 1.3 m × 0.7 m × 1.2 m. After the excavations, archaeologists identified it as an alabaster statue of Queen Tiye, the wife of Amenhotep III. Recently, Porcelli et al. [8] combined data from three geophysical methods (ERT, GPR, and magnetics) along with modern geomatic techniques to create an atlas of the Valley of the Kings, one of the most famous necropolises of Ancient Egypt. Ulrich et al. [9] applied the ERT method to an archaeological excavation at Tell Jenderes (northern Syria) and complemented the results of earlier geomagnetic and GPR studies by detecting Bronze Age settlement structures at greater depths than those methods could reach.
The two active-source electrical methods, ERT and GPR, complement each other regarding the information they provide for interpretation [10]. In this study, which aimed to detect void spaces, ERT estimates the volumetric properties of materials by measuring electrical resistivity or conductivity based on the spatial distribution of electrical potentials. The contribution of GPR and UST is the precise identification of boundaries between air and stone due to their significant variations in dielectric permittivity and acoustic impedance, respectively. Of the two methods, UST is much more sensitive to voids in solid materials because the impedance contrast is four magnitudes, while the permittivity contrast is only one magnitude (or half). Furthermore, and providing that the densities are different, UST can differentiate between materials with similar electrical or dielectric properties, which can be challenging for GPR and ERT alone [[11], [12]].
In this study, the combination of the three NDT techniques, i.e., ERT, GPR, and UST, was advantageous due to the suboptimal conditions for applying each technique individually on the Eastern face of the Menkaure Pyramid. ERT's major challenge was to overcome the high contact resistance of electrodes due to the high resistivity of the Aswan granite used for the Pyramid's casing, which has values in the tens of thousands of Ω [13]. For GPR and UST studies, limitations were associated with the presence of numerous joints between stone blocks, which, as a source of reflections, cause significant signal losses and thus reduce the depth of investigation. In addition, the depth of investigation was influenced by the shaling of some granite blocks.
First, an ERT survey was conducted over a large section of the Eastern face (26.95 m × 1.50 m) to identify high resistivity anomalies that may indicate the presence of air-filled voids. Subsequently, GPR and UST measurements were focused on anomalies observed in the ERT images to validate and accurately determine the position, size, and depth of these anomalies. The GPR survey was applied on a focused area (6 m × 4.05 m) using a 50 mm measurement grid to provide high resolution. Additionally, the UST measurements were applied over an even more focused area (5.40 m × 1.30 m) to confirm the previous ERT and GPR findings. IF was then used to combine the ERT, GPR, and UST images into a single composite image, allowing for the join interpretation of all measurements. Furthermore, a 3D laser-scanned model was used to create a sketch of the Eastern face of Menkaure to correlate the measurement locations of all the techniques used.
2. Methods
2.1. Electrical Resistivity Tomography (ERT)
2.1.1. Theoretical background of the ERT method
The ERT method is based on measuring of electrical resistivity, the intrinsic material property that describes the ability to withstand the flow of electrical current. ERT indicates differences in the electric properties of materials, thus enabling the recognition of resistive or conductive structures in relation to the surrounding medium. Measurements with the most frequently used four-point probe require two pairs of electrodes: two current electrodes - to feed the electrical current with a controllable strength and two electrodes installed away from the power source - to measure the voltage differences. Since such a configuration uses different electrode pairs to supply current and measure potential, contact resistance between electrodes and the sample has negligible influence on the measurement results. The resistance calculated from the applied current and measured voltage values depends on electrode network geometry. The value of the geometric factor unique to each electrode arrangement quantitatively specifies this relationship. The apparent resistivity values calculated with the contribution of the geometric factor are independent of the measurement geometry and represent the basis for a reconstruction of true resistivity distribution within the survey area during the iterative inversion process in which the electric potential distribution in the model space is computed using finite element (FE) or finite-difference (FD) solvers. In the first step, the inversion algorithm calculates the forward response of the reference model, which relies on known information, or assumptions, about the survey area and information about the measurement configuration. In the next step, the inversion program compares the predicted data with the observed apparent resistivity data and iterates through different resistivity models until the differences between them are reduced to an acceptable level. Information about measurement errors is required to narrow down the solution space for the inversion problem since it defines the limits within which the inversion model should match the observed apparent resistivity data. More detailed information on the theoretical principles of the ERT method, including data acquisition, limitations, and inversion techniques, can be found in the literature [[14], [15], [16], [17]].
2.1.2. ERT data acquisition and processing
ERT measurements were carried out using the high-resolution ABEM Terrameter LS2 measuring device (Fig. 2a). For this ERT study, a modified version of flat electrodes, specifically thin stainless steel mesh electrodes, were used to supply current and to measure voltage (Fig. 2b). The light weight of the mesh electrodes and their flexibility enable their easy deployment on the inclined wall of the Pyramid. Contact resistances were reduced by increasing the electrode size to 0.20 m × 0.20 m and using a synthetic sponge moistened with dilute salt water in electrode locations. Measurements were performed with two electrode configurations, Dipole-Dipole (DD) and Multiple Gradient (MG). The data were measured along four horizontal quasi-parallel lines with a separation between lines of 0.5 m and an electrode separation of 0.75 m (Fig. 2c). The ERT lines were placed on two adjacent rows of facing blocks, two lines in each row. The data were collected with varying current strength between 0.1 mA and 200 mA.
Fig. 2. (a) ERT field data acquisition, (b) close-up view of used stainless-steel mesh electrodes, and (c) layout of ERT lines on the sketch of the Eastern face of the Menkaure Pyramid.
In order to make a precise localization of the anomalies, a 3D CAD model (Fig. 3) was created with the guidance of a 3D laser scanner model measured by the ScanPyramids team. This ERT CAD model included two main solids: an inner region solid (inversion domain) with a size of 29.63 m × 5.46 m × 6.1 m (Fig. 3a) and an outer region solid. The 3D solids were combined in FreeCAD software [18] and meshed using Gmsh Software [19] (Fig. 3b). The final mesh file used in the inversion process consisted of 22,203 nodes, 117,807 cells, and 241,893 boundaries. All external boundaries of the 3D model were set to Neumann-type boundary conditions, so there is no current flow at the boundary to air, and the remaining boundaries were set to mixed-type boundary conditions.
Fig. 3. ERT data 3D model design (a) The designed 3D CAD inner region with layout of ERT lines, and (b) the 3D inner model after meshing in Gmsh (viewed in Paraview).
An open-source library, pyGIMLi, was used for geophysical inversion and numerical modelling [20]. To obtain a model of the Pyramid's face, a robust inverse solution L1 scheme with low regularization (λ = 5) and isotropic smoothing in different spatial directions was chosen. Compared to the L2 norm, which minimizes the sum of squared differences between observed and calculated data in a resistivity model, the L1 norm minimizes the sum of absolute differences, thereby allowing the selection of a model that is robust to outliers. The regularization strength (λ) is responsible for the model's complexity: an inversion with lower values of λ produces a model with more detail and hence sharper boundaries between different media [17]. The highest noise percentage was set to 4 % based on the reduced dataset's reciprocity measurement of the DD array. In addition, the absolute voltage-dependent noise was added to the data. Since lower voltages are generally measured with the DD array, the absolute value of the noise was chosen to be 20 μV. The root-mean-square (RMS) error was used to evaluate the differences between measured and predicted data in the inversion model.
The first step of the data processing was the data quality analysis procedure, which included analysis of stacking and reciprocity errors and checking for incorrect measured values. Due to time constraints during the surveys, reciprocal measurements were only collected for a reduced number of DD array datasets. The data showed varying reciprocal errors across different lines, ranging from 0.42 % to 5.45 %, with a mean value for all four lines of 3.55 %. Before inversion, the measured data were filtered to remove negative apparent resistivity values, as well as data with a stacking error exceeding 10 % and noticeable outliers. The 2D datasets were combined into a single 3D dataset and inverted using the pyGIMLi three-dimensional inversion algorithm.
Contact resistance checks were performed immediately before the start of each measurement. Most contact resistance values obtained from electrode tests were in the range from 1 kΩ to 15 kΩ, which is acceptable for high resistivity environments. They are close to values obtained under similar conditions, such as those reported by Pavoni, Carrera, and Boaga [21] in ERT measurements on wet rocky surfaces. The average contact resistance was 13 kΩ, and the highest contact resistance was 25 kΩ. The current strength in the survey ranged from 7 mA up to 130 mA, with an average value of 48 mA.
2.2. Ground Penetrating Radar
2.2.1. Theoretical background of the GPR method
GPR operates in the microwave range of the electromagnetic wave spectrum, which ranges from a few megahertz to 10 GHz. The method is based primarily on assessing the propagation of these waves through various host media. As soon as these waves encounter any heterogeneities (e.g., other media) or discontinuities (e.g., joint sets, faults, or subsurface utilities and voids), the waves are reflected and reach the surface in a few nanoseconds due to the high velocity. The reflected waves are collected by receivers in the acquisition equipment, and resulting radargrams are processed and interpreted to produce images of subsurface anomalies. First interpretations without post-processing can usually be made in near-realtime in situ. Anomalies could be discontinuities, changes in the rock properties, hidden voids and fractures, or groundwater tables.
Depending on the dielectric permittivity, which refers to the medium's ability to store a charge, waves travel through different media with different velocities. The permittivity of various mediums ranges from 1 (permittivity of air) to 81 (permittivity of water) and up to almost infinite values for metals. For the stone objects under investigation, permittivity values could generally range between 6 and 10, depending mainly on rock type and moisture content. The permittivity and the propagation velocity of the electromagnetic wave are inversely related: the higher the dielectric constant, the slower the propagation velocity in the medium. As a result, GPR reflections result based on the dielectric contrast between the host and reflector media, that is, the greater the contrast, the stronger the reflections. The recorded GPR data provides information about the types of anomalies and their depths. GPR applications span a wide range of disciplines and a lot of information can be found in the literature [[22], [23], [24], [25], [26], [27], [28], [29], [30]].
2.2.2. GPR data acquisition and processing
During the measurements, the zone where the ERT technique detected potential anomalies behind the granite casing of the Pyramid's Eastern face was scanned with dense profiles. The objective was to confirm these anomalies' existence and precisely determine their dimension. A RIS HI-MoD GPR system from IDS equipped with a 200/600 MHz dual frequency antenna was used. To perform the measurements in the most effective way, a customized cart was used along with a wooden bar to ensure the horizontality and verticality of the lines and to ease antenna movement on the face (Fig. 4a). Additionally, a wooden scaffold was constructed on the Eastern face of Menkaure to facilitate the measurements in the area of interest (Fig. 4b). A 50 mm grid was constructed to ensure sufficient resolution allowing for 3D image reconstruction (Fig. 4c).
Fig. 4. (a) GPR field data acquisition using 200/600 MHz by IDS system, (b) the built wooden scaffold in the area of study (c) layout of the focused 50 mm GPR grid measurements on the sketch of the Eastern face of the Menkaure Pyramid.
3D image reconstruction of the GPR datasets was carried out aiming to determine the shape and dimensions of these anomalies, which in turn improved the visualization. To construct 2D amplitude time slices, an interpolation scheme with linear weight for freely distributed 2D lines was adopted. The interpolation step was 50 mm in both the X and Y directions, and the trace increment, which mainly depends on the data acquisition system, was set to 13.2 mm. The time window was truncated at 50 ns because the signal had attenuated significantly after this value. The propagation velocity of the EM waves obtained using the depth-to-known reflector approach [32] was 0.13 m/ns (Fig. 5), which is within the range of expected values for granite found in the literature [33].
Fig. 5. An example of GPR Velocity Calibration on site, (a) a fallen granite stone at the site with a thickness of 0.76 m, (b) the corresponding GPR profile revealing a calculated velocity of 0.13 m/ns (after time conversion).
2.3. Ultrasonic testing
2.3.1. Theoretical background of the UST method
Ultrasonic testing (UST) is typically used to evaluate the internal structure and integrity of materials. The method has become an invaluable tool in the field of archaeology for non-invasive investigation and analysis of artifacts, structures, and monumental sites. UST operates within the frequency range of the acoustic wave spectrum, typically between 20 kHz and 200 kHz when applied to stone materials. These waves propagate through various materials at velocities depending on the material's density and their elastic properties.
As ultrasound waves travel through a medium, they encounter different features such as porosity, material interfaces, or internal defects. When the waves encounter such heterogeneities (e.g., in the form of different material layers) or discontinuities (e.g., cracks, voids, or inclusions), part of the wave energy is scattered or reflected back to the surface. An array of piezo-electric receivers collects the reflected waves. The data are then processed and reconstructed. Several methods are employed to reconstruct images from the processed ultrasonic data, for example, B-, C-Scan, Synthetic Aperture Focusing Technique (SAFT) [34], and Total Focusing Method (TFM) [35], which give the possibility of additional data processing changing the opening angle (aperture), amplitude (time gain), etc. Reconstructed images are processed and interpreted to deliver precise and dependable visual representations of the internal features of artifacts or subsurface structures.
The advantages of UST in archaeology using sensor arrays include its non-destructive nature, the ability to penetrate relatively deep into materials (up to 3 m), and its high sensitivity to small air-filled defects. Much like GPR, UST relies on the contrast between different material properties to detect anomalies. A comprehensive discussion of UST theory of sensing techniques, components, and instrumentations, as well as principles of data acquisition and processing, can be found in the literature, e.g., [34,36,37].
2.3.2. UST data acquisition and processing
UST measurements were conducted on the Eastern face of the Menkaure Pyramid with the main goal of confirming the anomalies detected by the ERT and GPR. To deal with the notable surface roughness, UST measurements were conducted using a sensor array with spring-loaded dry point-contact transducers that can handle a surface roughness of up to 6 mm without a need for a coupling agent. The UST field measurements were conducted using two 8-channel ultrasonic pulse-echo units (PD8000, Screening Eagle Technologies), combined to form a 16-channel array (Fig. 6a) at a total number of four horizontal profiles and 11 vertical profiles (Fig. 6b). The reflective scales mounted to the surface seen in Fig. 6a in combination with a built-in camera, allow for an automatic determination of the instrument's location with reference to the local coordinate system.
Fig. 6. (a) UST field data acquisition using two combined PD8000 arrays (16 channels), (b) the layout of the UST measurements on the sketch of the Eastern face of the Menkaure Pyramid.
The PD8000 has 24 dry point-contact (DPC) shear-wave transducers (three transducers in a row forming one channel) that act as emitters and receivers. One-by-one and left-to-right, each channel emits an ultrasonic pulse into the material while all subsequent channels record the response (Fig. 7). Using two combined units provides higher data coverage and penetration depth but also increases the weight of the instrument to be handled, stabilization time, and higher requirements of maintaining stable contact of all transducers with the surface. However, the number of individual signals for a single measurement increases from 28 to 120 while the recording time stays basically the same, which is a few seconds.
Fig. 7. Measurement process of PD8000 illustrated for two units combined (16-channel array).
For the measurements described in this article, a pulse frequency of 25 kHz and a waveform recording time of 4000 μs were chosen to achieve maximum penetration depth, which was found to be approximately 3 m. The measurement increment, i.e., the distance between two measurement locations in the X-direction, was defined as 0.10 m, and the grid spacing in the Y-direction was 0.15 m. To produce focused images, the shear (s-) wave velocity, which can be determined by using the instrument's calibration function, must be provided. A number of individual measurements at various locations of the inspection area were collected for this purpose. An average value of the shear wave velocity for the trapezoid block of 3400 m/s and for neighbouring blocks of 2815 m/s was determined from the measurements and stored for further data processing. For the measurements carried out in this research, reconstruction of the 2D images was achieved by further processing the recorded A-Scans via the Pundit Vision Software, which uses the Synthetic Aperture Focusing Technique (SAFT) [38].
3. Results
3.1. Electrical Resistivity Tomography (ERT)
After processing and 3D inversion of the ERT data set, the final true resistivity model was obtained and prepared for interpretation. In particular, to facilitate the analysis of the results, several 2D cross-sections in different directions were obtained from the inverted 3D model (Fig. 8a).
Fig. 8. (a) Relative location of the chosen 2D X and Y sections, (b) inversion results of the Y section with the presence of anomalies A1 and A2, (c) inversion results of the X and Z sections with the presence of anomalies A1 and A2, and (d) anomalies A1 and A2 overlayed on the Eastern face Menkaure sketch.
As can be seen from Fig. 8b (in section Y1), two layers with different background resistivities can be distinguished in the section. The upper layer with a thickness of 1 m–1.5 m has high resistivity values in the range of 5000 Ωm to 20,000 Ωm, and the underlying more conductive layer has resistivity values in the range of 500 Ωm to 1000 Ωm. Within the surface layer, the resistivity values are non-uniform. Thus, the inverted resistivity model generally correctly depicts the structure of the Pyramid face, allowing the separation of a high-resistivity (5000 Ωm to 20,000 Ωm) layer of granite blocks and the underlying layer of limestone blocks (500 Ωm to 1000 Ωm).
In addition to that, on sections Y1, X2, X1, Z1, and Z4 (Fig. 8b–c), two zones of anomalously high resistivity are visible, both of which lie behind blocks of polished granite. The first anomaly (A1) (>35,000 Ωm) is located behind the trapezoidal block, starting directly from the surface and located at a distance between 0.5 m and 1.8 m along the x-axis, while the second anomaly (A2) starts at 0.8 m depth and is located between x = −0.8 m and −2.8 m. Both anomalies are still visible at a depth of 2 m (Fig. 8c) but become less intense and smaller in size (Fig. 8c). The relative location and sizes of the anomalies on the Eastern face's sketch are presented in Fig. 8d.
Several possibilities should be considered to explain the high resistivity anomalies: 1) granite with physical properties that differ from neighbouring granite blocks (i.e., less porous and denser), 2) large air-filled gaps at the boundary between stone blocks, and 3) air-filled void. The most likely reason for both anomalies appears to be either an air-filled void or a large gap between the granite blocks of the outer casing and the limestone blocks.
3.2. Ground Penetrating Radar
After processing and 3D construction of the GPR data, the final GPR results were obtained and prepared for interpretation. In this section, the results from the 200 MHz frequency antenna are presented where three GPR profiles in the two measuring directions were displayed from the measured 3D GPR grid (Fig. 9). The results of the 600 MHz frequency are presented in Appendix A.
Fig. 9. (a) Layout of GPR profiles with marking 3 chosen profiles, (b) X direction 2D profile S1, andY direction 2D profile S2 and S3, (c) Time slices showing the relative location of A1 and A2 overlayed on the Eastern face Menkaure sketch.
The processed GPR images confirm the presence of the two anomalies (A1 and A2), which were detected previously by the ERT method, adding up that both anomalies are slightly inclined (Fig. 9). The first anomaly (A1) is located in the center of the area measured by GPR, behind the trapezoidal granite block, and the second anomaly (A2) is in the upper left part of the Pyramid face. The two anomalies differ in intensity and start at different depths.
Three radargrams were selected to display the GPR results with the detected anomalies (Fig. 9a). Profile (S1), in the X direction, revealed the presence of two inclined anomalies, A1 and A2 (Fig. 9b), while the other two sections (S2, S3) are oriented perpendicularly, intersecting A2 and A1, respectively (Fig. 9b). The time slice analysis of the results shows that the inclined anomaly (A1) began at a depth of 1.4 m, while the starting depth of anomaly (A2) was 1.14 m (Fig. 9c). Time slices were superimposed together to determine the average sizes of anomalies (Fig. 9c). They can be estimated as 1.5 m × 1 m and 0.9 m × 0.7 m for A1 and A2, respectively.
3.3. Ultrasonic Testing
After processing UST data and reconstructing images, various reflections were observed. Two main anomalies were localized by the UST measurements, confirming the detected anomalies. One dominant reflector lies behind the trapezoid granite block (A1), in addition to the presence of a second reflector at the same location as detected anomaly A2. The two detected anomalies exhibit different attributes concerning their depth and extension. Fig. 10 shows a reconstructed image of the horizontal profile H02 (Fig. 10a–b), which cuts through the two anomalies and vertical profiles V02 and V06, defining the variable amplitude of the beginning and end position of anomaly A1 (Fig. 10b–d).
Fig. 10. (a) The location of the profiles of the reconstructed UST image (in blue), (b) the reconstructed UST image of H02 profile, (c) the reconstructed UST image of V02 profile, and (d) The reconstructed UST image of V06.
The dominant reflector with a hyperbolic-like inclined shape can be observed over X = 2.6 m–4.1 m (A1). This type of response can be caused by a change in the material or a discontinuity on the surface of the backwall. Considering data from vertical and horizontal profiles, it can be assumed that anomaly A1 is located at a depth between 1.3 m and 1.6 m.
The second reflector with a hyperbolic-like inclined shape, observed across X = 0.64–1.9 m, matches anomaly A2 and is detected at an average depth of 1.12 m. Due to the roughness of the surface and the presence of near-surface shales, measuring vertical profiles in this block was not possible. On the other hand, the multiple inclined reflectors (J), starting from the surface and leading through the entire profile, are caused by the presence of block joints in the horizontal direction.
3.4. Image fusion
Image fusion (IF) was employed to merge the images discussed in the previous sections and interpret them. The general procedure used herein is presented in Ref. [5] is described briefly step-by-step in Appendix (B). IF was performed on two selected profiles illustrated in Fig. 11a.
Fig. 11. (a) The sketch of studied area showing profiles P1 and P2 and coordination origin, (b) final fused image for horizontal profile P1 (Y = −0.20 m) with the red line marking profile P2 shown in Fig. c, (c) final fused image for vertical profile P2 (X = +1.40 m) with the red line marking profile P1 shown in Fig. b. Note that the extent of the ERT image is limited to Y ≈ −1.00 m. The input and fused images can be found in Appendix (B), Figs. B.17 and B.18.
Fig. 11b shows the final fused image for the horizontal profile P1 at Y = −0.20 m. Two strong inclined reflectors are visible in the image, labelled 'a1' and 'a2'. While reflector a1 lies behind the high resistivity region of anomaly A1, reflector a2 is shifted in the X-direction relative to anomaly A2. Note that anomaly 'A1' can be explained by the trapezoidal block's (highlighted in light blue) high resistivity compared to the neighbouring blocks, which is discussed in Sections 3.1 Electrical Resistivity Tomography (ERT), 4 Discussion. Due to its deeper location, Anomaly 'A2' can be associated more likely with an air-filled void. In Fig. 11c, the final fused image for the vertical profile P2 at X = 1.40 m is shown. The strong diagonal reflector labeled 'b' is visible in both the GPR and UST images and can be associated with the joint between the trapezoidal block and the block supporting it. Reflector 'a1' corresponds to reflector 'a1' in Fig. 11b. A fused image has several advantages over any individual input image. First, the variables that must be selected to create the input images for IF (e.g., shear wave velocity for UST image) are confirmed when the reflectors appear in the same position in the fused image, which is the case here. Second, each technique has its strengths and limitations and combining them helps overcome the limitations of each individual technique. As an example, the unique strength of ERT is that it shows a volume rather than an internal boundary, which, on the other hand, GPR and UST images provide. Finally, pertinent information from all input images are contained in a single composite image, which is easier to interpret by a human inspector as well as an algorithm. The GPR, UST and ERT input images used to create the final fused images for Profiles P1 and P2 are shown in Appendix (B), Fig. B.17 and B.18, respectively.
4. Discussion
Over the past three years, successive measurements using ERT, GPR, and UST techniques have been conducted to identify potential anomalies behind the Eastern face of the Menkaure Pyramid. These measurements aimed to test the hypothesis of a second entrance (void) behind the granite casing blocks.
In the initial stage, the ERT method identified two high resistivity anomalies, A1 and A2, from a large measured area of the Eastern face, measuring 26.95 m × 1.5 m. The next stage involved a detailed study using GPR and UST methods at the predetermined positions based on the ERT results. In contrast to ERT, the GPR and UST methods were able to accurately determine both the starting depths and the dimensions of the anomalies. Furthermore, they also reveal that the anomalies are slightly inclined. Anomaly A1 measures 1.5 m × 1.0 m at a starting depth of 1.4 m, and anomaly A2 measures 0.9 m × 0.7 m at a starting depth of 1.13 m. These starting depths were calculated as the average of the GPR and UST data, while the x and y dimensions were based on GPR time slice results.
The final stage of the study was to conduct IF to compare the results of all three methods and, more clearly, to highlight the similarities and differences between the detected anomalies by each technique regarding position, shape, inclination, and depth. As shown by the image fusion results, the detected position of the anomalies from all three methods corresponded (Fig. 11). However, there were slight differences in the depths at which the anomalies appeared, along with a minor shift in the X direction of anomaly A2. According to GPR and UST data, anomaly A1 begins at a depth of 1.4 m, and anomaly A2 begins at a depth of 1.13 m. In contrast, the ERT data show that anomaly A1 starts directly from the surface and extends to a depth of at least 2 m, while anomaly A2 begins at a depth of approximately 0.80 m and has a similar extent. The fact that anomaly A1 started directly from the surface in the ERT results can be explained by the high resistivity of the trapezoidal block, which was detected by the relatively high contact resistance of the trapezoidal block compared to the surrounding blocks during the electrode test. In such a case, the contrast between the air-filled void and the high-resistivity granite block may not be sufficient for a clear separation of anomalies in the ERT inversion results. Since anomaly A2 is located at the edge of the GPR grid, it is recommended to consider both the GPR/UST and ERT locations for a joint interpretation and localization.
Various factors that represented challenges for interpreting the ERT data influenced the electrical resistivity distribution in the study area. These factors included the presence of several types of material (limestone, granite, and filling material), joints and air gaps between stone blocks, and variations in weathering degree between the superficial and deeper stone layers. Additionally, the presence of the high resistivity granite limited the total depth of investigation. The main difficulty in interpreting GPR data—determining the extension of the anomalies — is caused by the rapid decrease in signal strength with depth. Similarly, UST was affected by the high attenuation of the granite stones and their deterioration. This, along with the rough surface, made it not only difficult to determine the shear wave velocity but also to carry out the measurements.
Table 1. The physical parameters of the materials involved in the simulations.
Material
Resistivity (Ω m)
Dielectric Const.
Shear Velocity (m/s)
Density (kg/cc)
Granite
20,000
5
2815
2750
Resistive granite block (Trapezoid block)
40,000
4.5
3400
3000
Limestone
1000
7.5
2400
2711
Air-voids
1,000,000
1
0
1.2
From the simulation results, it was found that if there is a long air-filled void behind a high resistive granite block (higher than the surrounding medium), the boundary between both mediums cannot be resolved by ERT while it can be clearly detected using GPR and UST. Furthermore, it was clear that anomaly A2 fits the possible scenario of the presence of a simple air-filled void (Model 1, Appendix C), but it is not clear what is the exact length of such a void. For anomaly A1, further analysis was required due to the complexity of the ERT field results.
After carrying out more extensive 2D ERT, GPR, and UST simulations and analysing their results, it was found that the most plausible scenario that could explain anomaly A1 is the presence of a layer of granite with one high resistive granite block (the trapezoid block) followed by an air-filled void within the limestone layer (Fig. 12a). The simulation results are shown in Fig. 12b–d and are comparable to the field results of ERT, GPR and UST.
Fig. 12. (a) Sketch of the designed simulation, (b) Results of ERT simulation, (c) Results of GPR simulation, (d) Results of UST simulation.
Despite the adopted approach successfully detecting two anomalies (Fig. 13), the extension of these anomalies is yet to be determined. Consequently, further investigations are necessary using ERT, GPR, and UST simulations and physical models. In addition to include other techniques such as infrared thermography (IR), microgravimetry, and cosmic-ray muons to give more information regarding the extension and the possible explanations for these anomalies.
Fig. 13. The location and dimensions of the detected anomalies, (a) on the sketch of the Eastern face of Menkaure and (b) on an actual photograph of the Eastern face of Menkaure.
5. Conclusion
An integrated approach combining three NDT techniques (ERT, GPR, and UST) was utilized along with multimodal IF. The goal was to investigate the Eastern face of the Menkaure Pyramid, with a particular focus on identifying anomalies behind the polished granite blocks. The interpretation was supported by 2D numerical simulations of GPR, UST, and ERT for selecting the possible scenarios that explain the actual measurements. Furthermore, this integrated approach resulted in overcoming the limitations regarding the application of a single technique and enhancing the validity of the interpretations for the complex medium.
After processing and applying the joint interpretation to the field data, two main anomalies, A1 and A2, were identified behind the granite blocks of the Eastern face. It was found that anomaly A1 has the approximate dimensions of 1.5 m × 1.0 m and starts at a depth of 1.35 m, while the approximate dimensions of anomaly A2 are 0.9 m × 0.7 m starting at a depth of 1.13 m. Based on extensive 2D ERT, GPR, and UST simulations, as well as IF, it was concluded that both anomalies likely represent air-filled voids within a limestone medium, beginning directly behind the outer granite casing of the Eastern face. Notably, anomaly A1 is located behind a high resistive trapezoid granite block.
Although the integrated approach was successfully able to detect the approximate dimensions and the starting depths of the two potential anomalies, it was difficult to determine how far the anomalies extend inside the pyramid due to limitations in the penetration depth of the used methods.
It becomes clear that each technique used in this study has its own limitations. Consequently, the combination of the three techniques gave a better understating of the detected anomalies. However, the interpretation of the detected anomalies should be discussed among Egyptologists. Additionally, we recommend further investigation of the anomalies using other non-destructive techniques to overcome the encountered limitations of the study.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
We acknowledge the Science, Technology & Innovation Funding Authority (STDF) for funding the establishment of a new center of excellence in Non-destructive Techniques & Engineering Geophysics at Cairo University - Faculty of Engineering - (project ID 43842) providing the equipment used in this research. The authors acknowledge the support of yTUM International Graduate School of Science and Engineering (IGSSE) as well as the German Academic Exchange Service (DAAD) in the scope of the "German-Egyptian Progress Partnership, Program Line 2″ under the title "Non-Destructive Techniques for the Preservation of Egyptian Cultural Heritage". Special thanks are due to the Supreme Council of Antiquities and the Egyptian Ministry of Tourism and Antiquities for their support of the ScanPyramids mission. We thank Mr. Sherif Fathy (Minister of Tourism and Antiquities), Mr. Ahmed Issa, Dr. Mamdouh Eldamaty and Dr. Khaled El-Enany (Former Ministers of Tourism and Antiquities), Dr. Mohamed Ismael (Secretary General of the Supreme Council of Antiquities), the Scientific Archaeological Committee, Mr. Ashraf Mohy (Director of the Pyramids Archaeological Area), and their collaborators and assistants. Thanks to the management of Cairo University and the Faculty of Engineering - Cairo University for the support of the ScanPyramids project. Special thanks to Mr. Johannes Scherr (Technical University of Munich) for his support in the first measurements campaign and Hamada Anwar, the logistics coordinator of site and international missions of the ScanPyramids project. We would like to sincerely thank Dr. Norbert Klitzsch for his consistent support throughout the ERT study. We are also grateful to Dr. Carsten Rücker and Prof. Dr. Thomas Günther for their helpful advice on working with the pyGIMLi software. The result of this research is part of the ScanPyramids project, which is supported by: NHK, la Fondation Dassault Systèmes, The French Embassy in Egypt, TNG Technology Consulting, and Mondaic AG, which have provided unfailing support to the project.
Appendix A. Additional Ground Penetrating Radar Results:
A 3D analysis of the GPR datasets was carried out to determine the shape and dimensions of these anomalies and to improve visualization. To construct a three-dimensional image of an iso-amplitude time slice, an interpolation scheme with linear weight for freely distributed 2D lines was adopted. The interpolation step was 50 mm in both the X and Y directions, and the trace-increment, which mainly depends on the data acquisition system, was set to 13.2 mm. The time window was truncated at 50 ns due to significant signal attenuation beyond this point. A summary of the GPR results of the two frequency channels is presented in Figs. A.14–15 showing the time slices of for A1 and A2 respectively and the corresponding depths.
Fig. A.14. Time slices and their corresponding depths of anomaly A1 overlayed on the Eastern Face laser scanner sketch (a) for 200 MHz. channel, (b) for 600 MHz. channel.
Fig. A.15. Time slices and their corresponding depths of anomaly A2 overlayed on the Eastern Face laser scanner sketch (a) for 200 MHz. channel, (b) for 600 MHz. channel.
Appendix B. Additional Information about Image Fusion:
This section discusses the image fusion (IF) procedure used in this study (see Fig. B.16) and shows all input and fused images for Profiles P1 and P2 (see Figs. B.17 and B.18, respectively). Pixel-level IF using a Discrete Wavelet Transform (DWT)-based algorithm was used to merge the reconstructed GPR, UST, and ERT images into a single composite image. The GPR and UST images (see Figs. B.17/B.18a and b, respectively) were created using synthetic aperture focusing technique (SAFT) and total focusing method (TFM)-based algorithms (see Sections 2.2.2 GPR data acquisition and processing, 2.3.1 Theoretical background of the UST method, respectively) and show the reflections from internal boundaries between different materials. The ERT image (see Figs. B.17/B.18d) is a tomographic slice that provides information regarding the resistivity distribution (see Sections 2.1.1 Theoretical background of the ERT method, 2.1.2 ERT data acquisition and processing). The GPR and UST images were kept in grayscale and the ERT image was defined using the "jet" colormap. Note that only two images were fused at a time, requiring a two-step IF process. First, the GPR and UST images were fused, creating the "Fused GPR/UST" image (see Figs. B.17/B.18c). Subsequently, this image was merged with the ERT image, which is a tomographic resistivity slice, to produce the "Final fused" image (see Figs. B.17/B.18e). Following Fig. B.16, the IF procedure shall be briefly described next.
The following steps were applied to the individual input (i.e., GPR, UST, and ERT) images:
•
Registration – This step involves translating and cropping the input images so that they are properly aligned and contain the same content, making sure the images are compatible.
•
Resampling – In order to perform pixel-level IF, the input images must have the same dimensions (or same number of pixels in the x and y directions), which was achieved by resampling them using the MATLAB function imresize().The image resolution should be selected so that the smallest features (or wavelengths) of any given input image is properly sampled. In this study, a pixel size of 2 mm × 2 mm was selected as it provided sufficient resolution to capture the shortest wavelengths which are present in the UST image.
•
Intensity scaling – This step was used to normalize the intensity values so that they take values from 0 to 1, meaning that each input image has the same weight (or importance). Prior to normalization, MATLAB's imadjust() function was used to increase the contrast of the GPR and UST images. The intensity of the ERT image did not require contrast adjustment.
Fusion - DWT-based IF was employed using MATLAB's wfusimg () function [43,44], which decomposes the two input images using specified fusion rules into approximations and details coefficients. The "sym3" wavelet with a wavelet decomposition level of seven was selected for all processing and fusion was performed by taking the minimum and maximum for the approximations and details coefficients, respectively. These settings were determined iteratively and manually (as illustrated in Fig. B.16) until found to produce images with optimal contrast.
The final step is the interpretation of the final fused images (see Figs. B.17/B.18e), which was done by working collaboratively as a multi-disciplinary team.
Figs. B.17 and 18 show the input as well as fused images for Profile P1 and P2, respectively (locations see Fig. 11a). Note the clearly visible reflectors in the GPR and UST images (Figs. B.17/B.18a and b, respectively) as well as the high resistivity regions in Figs. B.17/B.18d. The fused GPR/UST images (Figs. B.17/B.18c) contain pertinent features from both input images with the two strong reflectors being aligned exactly. The final fused images (Figs. B.17/B.18e) allow a holistic view of important reflectors as well as regions of different resistivity. A detailed interpretation of Figs. B.17e and B.18e is provided in Section 3.4 for Fig. 11b and c, respectively.
Fig. B.17. Input and fused images for final fused image shown in Fig. 11b (Profile P1 @ Y = −0.20 m): (a) GPR (200 MHz), (b) UST, (c) fused GPR and UST, (d) ERT, and (e) final fused image (see Fig. 11b). IF = image fusion.
Fig. B.18. Input and fused images for final fused image shown in Fig. 11c (Profile P2 @ X = 1.40 m): (a) GPR (200 MHz), (b) UST, (c) fused GPR and UST, (d) ERT, and (e) final fused image (see Fig. 11c). IF = image fusion.
Appendix C. Forward modelling of ERT, GPR, and UST methods:
of 15 and the isotropic smoothness constraint (zWeight = 1) were used to calculate all forward models. A very fine mesh was chosen for the forward models with a maximum triangle mesh size of 0.01 m2.
The ERT simulations were conducted along 26.86 m 2D profiles while GPR and UST simulations were conducted along 6 m 2D profiles (similar dimensions and locations for the field-measured profiles). The main materials involved in the simulation and their physical parameters are illustrated in Table 1.
Five main scenarios were tested and numerically simulated to study the possible interpretation of the measured field data, and they are as follows:
The first simulation was carried out for a simple case with the first layer of granite, followed at a depth of 1.35–1.61 m by a limestone layer (Fig. C.19). To investigate the possibility of the existence of a potential corridor, an inclined air-filled void was simulated at a depth of 1.4–1.55 m, started directly after the end of the trapezoidal granite block, and continued to 3.5 m depth. This void was started by x = 0.7 m and ended by x = 1.7 m. The second scenario simulated several 5 mm joints between the granite blocks followed by a homogenous limestone layer (Fig. C.20). Scenario 3 exhibited a high resistive trapezoid block (denser with higher seismic velocity) inside the granite layer followed by a homogenous limestone layer (Fig. C.21). In scenario 4, the high resistive trapezoid block, introduced in the scenario 3, along with several 5 mm joints between the granite blocks followed by a homogenous limestone layer was presented (Fig. C.22). The fifth scenario was illustrated in the main text (Fig. 12).
Fig. C.19. (a) Sketch of the designed simulation (scenario 1), (b) Results of ERT simulation, (c) Results of GPR simulation, (d) Results of UST simulation.
Fig. C.20. (a) Sketch of the designed simulation (scenario 2), (b) Results of ERT simulation, (c) Results of GPR simulation, (d) Results of UST simulation.
Fig. C.21. (a) Sketch of the designed simulation (scenario 3), (b) Results of ERT simulation, (c) Results of GPR simulation, (d) Results of UST simulation.
Fig. C.22. (a) Sketch of the designed simulation (scenario 4), (b) Results of ERT simulation, (c) Results of GPR simulation, (d) Results of UST simulation.
Data availability
Data will be made available on request and after permission from the Egyptian Ministry of Tourism and Antiquities.
Shallow seismic refraction, two-dimensional electrical resistivity imaging, and ground penetrating radar for imaging the ancient monuments at the western shore of Old Luxor city, Egypt
B. Riveiro, M. Solla (Eds.), Non-destructive techniques for the evaluation of structures and infrastructure, vol. 11, CRC Press, Boca Raton, FL, USA (2016)
No comments:
Post a Comment