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)
