Integrated Approach to Investigating Historic Cemeteries

: Ground-penetrating radar (GPR) and terrestrial laser scanning (TLS) surveys were conducted at a historic cemetery at Cape Canaveral Air Force Station, Florida, U.S., in order to conﬁrm the presence of burials corresponding to grave markers and detect potential unmarked burials. Noise in the GPR data from surface features and subtle terrain differences must be addressed to determine the extent of anomalies of interest. We use singular value decomposition (SVD) to isolate and remove energy from GPR data. SVD allows one to remove unwanted signals that traditional processing techniques cannot. With SVD ﬁltering, we resolve an anomaly adjacent to conﬁrmed burials otherwise overprinted by unwanted signal. The migration of SVD-ﬁltered data produces more distinct, spatially constrained point reﬂectors. Ground elevation is derived from georeferenced TLS data and compared to that from airborne laser scanning (ALS) to highlight subtle terrain that can assist data interpretation. TLS elevations show a subtle modern mound over the burial plot where ALS elevations show a depression. The targets of interest are approximately 20–30 cm higher in elevation if a topographic correction is performed using TLS versus ALS. In archaeological applications, a notable change is often recorded at the sub-meter scale. The combined approach presented here better resolves geophysical response of buried features and their positions in the ground relative to each other.


Introduction
Ground-penetrating radar (GPR) has long been used for determining the location, depth, and shape of human burials e.g., [1][2][3][4][5][6]. Because of its non-invasive approach, rapid data acquisition, and real-time data viewing, GPR is an ideal method for forensic and archaeological applications. Nevertheless, the method does suffer from several significant limitations, including, but not limited to, high electromagnetic signal attenuation (loss of signal with depth), scattering (energy reflected in multiple directions), and false positives (misidentification) from roots or other buried objects (e.g., [3,7,8]). Although GPR has been well studied and applied across many disciplines, extracting sharp images (highlighting desired anomalies while minimizing unwanted features) through signal processing remains a challenge.
The GPR response of a grave is sometimes a diffraction pattern (hyperbola-shaped signature often produced by point-like targets) or a disturbed section of strata. In both scenarios, interpretations can be limited. For example, migration algorithms collapse diffraction events to single points that are positioned at a depth corresponding to the diffraction apex to produce a more realistic representation of the subsurface [9]. However, other nearby buried features (i.e., cobbles, tree roots, or ringing from shallower features) are erroneously incorporated into the migration scheme. The result is anomalies that comprise more than one target, making it challenging to distinguish individual buried objects.
Singular value decomposition (SVD) is an advantageous filtering technique: reduce high-dimensional, highly variable data to a lower-dimensional space that reveals substructures, discard the unwanted substructures (i.e., direct wave, ringing, random noise), and reconstruct a filtered image. SVD has been previously applied to GPR data e.g., [10][11][12][13][14][15][16][17] with much success. In particular, SVD has been used to increase the resolution of small, point-like reflections in GPR [11]. The whole approach is computationally efficient and ideal for data sets in which the targets and non-targets have non-horizontal geometries [10] and contrasting amplitudes [17].
In areas with variable relief, a topographic correction of GPR data should be performed either as a topographic migration [18,19] or simply a static shift [20]. In areas of very low but observable relief, a static shift following migration is sufficient. Global Positioning Systems (GPS) can be paired with GPR equipment to synchronously collect geospatial and geophysical data. It has been shown that high-precision height information is valuable in archaeological settings e.g., [21][22][23][24][25]. Digital elevation models (DEM) derived from airborne laser scanning (ALS) datasets can provide full coverage of a survey area (ALS is also referred to as airborne Light Detection and Ranging or airborne LiDAR). Elevation can also be extracted from point cloud data e.g., [24][25][26][27] to perform a topographic correction of GPR data. Alternatively, elevation models derived from terrestrial laser scan (TLS) data provide a higher horizontal and vertical resolution that can reveal subtle mounds or depressions that may correspond to subsurface features [28].
This paper presents the geophysical findings from the Wilson Brothers Cemetery on Cape Canaveral Air Force Station, Florida (CCAFS, Figure 1a) and the joint data processing and remote sensing techniques used to increase GPR image resolution and refine the shape and depth of anomalies. We show that SVD can remove ringing that overprints targets of interest and improves migration results. We then perform topographic corrections while using elevation derived from TLS and ALS to highlight subtle topographic changes that can influence GPR data interpretation. The techniques themselves are well-known, but have not been jointly applied to detect and characterize burials.

Historical Context
Many Anglo-American homesteaders moved to Cape Canaveral after the Civil War to exploit expanded economic ventures in the fishing and agricultural industries. Eight cemeteries were established on the cape dating from the late 1800s to mid-1950s, mostly small family plots, one community cemetery, and two believed to be single burials. The United States government began purchasing land on the island in the late 1930s in order to establish a long-range proving ground for missiles. By 1952, the entirety of the cape was federal land. As a cape occupied for at least 3000 years [29,30], the remains of several distinct archaeological site types, including the eight historic cemeteries, are spatially intertwined with the iconic launch complexes and Air Force facilities that are embedded with historical significance.
Most of the cemeteries were left abandoned or poorly maintained until relatively recently. The chain-link fences that border each cemetery are not original and are considered arbitrary boundaries. The exact footprint of each cemetery is unknown, and many cemeteries were speculated to have unmarked burials. Therefore, in support of CCASF cultural resource management GPR, GPS, and TLS were undertaken at each cemetery in order to verify cemetery extent and locate any evidence of unmarked burials.
The Wilson Brothers Cemetery, recorded with the Florida Master Site File as 8BR3366, is a small Christian family plot that includes two known graves, Frank and Alfred Wilson. They are grandchildren of the locally well-known Mills Burnham, the first permanent lighthouse keeper on Cape Canaveral. Frank and Alfred were both born on the island and passed away in 1940 and 1942, respectively. Any specific religious denomination that would influence burial practices is unknown.
The fenced area is approximately 126 square meters. Within the fence, two east-facing headstones are enclosed by railroad ties that are approximately 15 cm × 15 cm × 2 m planks of wood ( Figure 1b). Shell midden material, attributable to an ancient archaeological site, is also present in this area, but appears scattered and ephemeral.
The railroad ties present a challenge not encountered at the other cemeteries. The GPR signal from the ties is high amplitude ringing-coherent noise where a reverberating response from an object oscillates far longer than the GPR pulse. That oscillation is recorded as repeating reflections when, in reality, it is the energy reflected off of a single object. This site required the use of filtering techniques that the remaining cemetery datasets did not.

Physical Context
The surficial geology of CCAFS comprises unconsolidated deposits of fine-to medium-grain sand, silt, shell, and clay over limestone [31][32][33]. The cemetery is located on the west coast of the cape along a north-south trending sand ridge. The environment is generally characterized by high infiltration rates [34]. Soil cores taken in the vicinity of the cemetery show varying ratios of shell, sand, and silt down to at least 7 m [35]. Beneath these sediments, clay content increases with depth. Bulk soil composition is an important consideration, because shells can produce numerous high amplitude reflections, and radar does not penetrate clay. The depth of investigation for this study does not intersect with clay layers. High shell content (natural and cultural) in the uppermost material is considered in data interpretations.

Materials and Methods
Common-offset (fixed distance between the transmitter and receiver) GPR data were collected with a 500 MHz shielded antenna and MALÅProEx system, where the system was mounted on a cart and pushed along the ground surface. Transects were collected in a grid and spaced 0.5 m apart and oriented perpendicular to presumed burial orientation. The trace interval (distance between scans) for all transects is 0.025 m. The cart odometer recorded transect distance. This level of spatial sampling does not obtain full resolution of the subsurface [36], but enough to resolve burial-sized targets and disturbed strata. Geophysical data are protected according to 2019 Florida Statutes Title XVIII, Chapter 267.
GPS information was collected from all surface features (i.e., trees, tree stumps, fence corners, grave markers) and the start and stop position of each GPR transect-defined as the center position of the antenna. Field collection was performed using a Trimble GeoExplorer GeoXH Network Rover with a Zephyr2 antenna compatible with the sub-meter and centimeter package post-processing.
A simple dewow and time-zero correction were performed on each profile before SVD to remove inherent low-frequency noise and designate the starting point of the received signal, respectively. Singular value decomposition was performed in MATLAB (see Appendix A for details). Reconstructed GPR data were transferred to GPR-SLICE for positioning, migration, three-dimensional (3D) volume interpolation, and topographic correction. We used hyperbola-fitting to estimate a single velocity of 0.16 m/ns, which reflect the dry conditions during data acquisition.
The 3D volumes were obtained in GPR-SLICE by sampling all profiles over the same increments of two-way travel time (2.9 ns) to produce 40 binned slices. The amplitudes of each bin are the average of values sampled within that bin. The slices were gridded at a 0.05 m cell size with inverse distance interpolation between samples of a given slice. The 3D volume is comprised of these grids and interpolations between them. Surface elevations that were derived from ALS and TLS were each imported for separate topographic corrections. Each dataset was gridded over the same cell size as the GPR data (0.05 by 0.05 m). The 3D volume was then draped over topography. Topographic profiles were also extracted to perform a static shift on select GPR profiles.
ALS data were published in 2007 by the Florida Division of Emergency Management [37]. These have 1 m horizontal accuracy and 0.18 m vertical accuracy. Only points classified as ground were included in producing the 30-cm DEM.
A hybrid phase-shift scanner (Leica RTC360) was used in the TLS survey. Eight scans were collected in the fall of 2019. Data capture included all surface locations within the cemetery boundary. The scans were registered using native to scanner Cyclone software. The registration error was less than 5 mm on average. The final point cloud has a minimum point spacing of 5 mm along the ground. A digital elevation model (DEM) was produced by georeferencing the TLS point cloud. The global horizontal positions of cemetery fence corner posts provide four reference points that were aligned with corresponding points in the TLS point cloud while using the open-source software CloudCompare. The corner elevations were extracted from the ALS-derived DEM. The wider scale terrain trends between ALS and TLS agree, so reference point elevations serve as elevation controls for otherwise arbitrary TLS height values. The main objective was to resolve small terrain differences within the cemetery. Once the TLS point cloud was georeferenced, all surface features (trees, fencing, headstones) were removed.
Both elevation datasets were gridded with 0.05 m cell size with inverse distance interpolation, which is the same cell size of the GPR depth slices. Because the ground is maintained grass, a minimum elevation value was assigned to each cell for the TLS data. The ALS-derived DEM has 30 cm resolution so the smaller cell size gained no resolution. TLS point clouds have subcentimeter scale resolution. Thus, the resolution was decreased by gridding to that of the GPR depth slices. Nevertheless, the detail captured by TLS and not ALS is still apparent.

Results
A total of 24 subparallel GPR transects (212.5 m) were collected over the Wilson Brothers Cemetery in a northeast-southwest direction. The closest soil core (approximately 1.5 km to the south) indicates 1.5 m of sand underlain by another 1.5 m of shell-rich sand. GPR data show the top of a semicontinuous layer of high amplitude reflectors at about 2.5 m depth, which is interpreted as the top of shell-rich sediments. Two anomalies that were identified as burials produce clear diffraction patterns in several successive profiles. Figure 2 presents an example profile crossing over the burials.
The partially buried railroad ties produced very strong ringing down to about a meter depth (Figure 2a, x = 4.5, 7 m). This ringing overprints a third diffraction pattern detected in several profiles. Although the diffraction is recognizable in profile view, it is difficult to resolve in depth slices due to the overprinting.

Singular Value Decomposition of GPR Data
We must remove the direct wave (the received pulse traveling directly to the receiver) and the ringing to reveal the depth and extent of the overprinted diffractions. A standard mean trace subtraction (a.k.a background removal) successfully removes much of the direct wave; however, the ringing is entirely untouched (Figure 2b). Rather than choosing a variety of trace widths through trial and error to calculate the best trace mean, we isolate the ringing and pull it out of the profile with SVD. A detailed explanation of the SVD approach is in Appendix A. SVD decomposes the GPR image into weighted eigenimages. To give the reader an idea of how much an image is decomposed, in the example profile in Figure 2, SVD returned 390 eigenimages. As expected, the direct wave is the most dominant structure and captured in the first eigenimage ( Figure A1a). The ringing produced by the railroad ties at the surface is captured in the second through fifth eigenimages ( Figure A1b). A reconstructed GPR image-summing all but the first five eigenimages-now reveals reflections coming from the bottom of the railroad ties and the third diffraction that may or may not be associated with the burial plot (Figure 2c). A large amount of energy is removed in the process, which is expected as the lowest eigenimage indices hold the most energy. Thus, an automatic gain control is applied after migration, but before creating a 3D volume and topographic correction.
In map view, ringing that is produced by the railroad ties is most apparent at the surface where it is essentially the only feature detected. (Figure 3a). SVD was performed, as described above, on all the GPR profiles before migration, gain, and topographic correction. The ringing is almost completely removed in map view (Figure 3b). At a depth of 0.5 m, the ringing is still apparent overprinting a shallow anomaly. Furthermore, the ringing and anomaly combined appear as a relatively wider section of the railroad ties (Figure 3c). With SVD performed on all profiles, the same shallow anomaly can be resolved as one or two features beneath the southern railroad tie (Figure 3d). A tree root system in the eastern half of the grid is also more apparent.

Migration of SVD-Filtered GPR Data
We performed a two-dimensional (2D) Kirchhoff migration before and after SVD using a constant ground velocity (0.16 m/ns calculated from diffraction hyperbola fitting). It has been shown that TLS can be utilized to determine a best fit ground velocity [24]. However, that is only possible when TLS data have captured subsurface structures (e.g., a cave or tunnel). We also acknowledge that a velocity profile would be more appropriate for resolving targets over a range of depths [19]. The depth of investigation in this study is only the first two meters, and a single velocity fits the diffraction patterns within that range.
Migration that includes the ringing has multiple issues. First, the velocity of the ringing differs from the ground where target diffractions are. The chosen velocity may appropriately collapse diffractions (Figure 2d, red circles), but it is an overestimate for the ringing. As a result, the signal of the ringing is spread out and obscures nearby features even more. Second, the ringing and diffraction pattern it overprints are collapsed together to an erroneously large single point reflector.
The migration of the SVD-filtered data produces more spatially constrained reflectors in the absence of ringing (Figure 2e). One can see what is interpreted as the bottom of the railroad ties at about 0.5 m below the surface. Additionally, a third point-like anomaly (Figure 2e, yellow circle) beneath the ties are observed where only two are discernible in the original data. The ability to migrate the two targets of interest at depth does not appear to be negatively affected by SVD.

Topographic Correction Using Elevations from ALS and TLS
Elevation comes from ALS-and TLS-derived elevation models (Figure 4a,b, respectively). The ALS DEM measures a depression east of the graves that is not detected by TLS. Rather, TLS shows a square-shaped mound corresponding to the graves' location, which is built-up gravel bounded by the railroad ties.  The variations between DEMs are small but enough to cause slightly different anomalies to appear in a given elevation slice. We illustrate this by presenting a pair of elevation slices for a target of interest where one dataset is topographically corrected with ALS-derived elevations and the other with TLS-derived elevations ( Figure 5). The targets are all approximately 20 cm higher in elevation when shifted by TLS values versus ALS values (Figure 5e,f). The exception is the northernmost burial (Figure 5g,h). This feature is 30 cm higher in TLS elevation slices, which is expected. This area corresponds to the edge of a depression in the ALS DEM, but not in the TLS DEM.
Beyond the targets of interest, there is disagreement between elevation slices. Tree roots extending into the survey area from the east and southeast appear to differing degrees between slices. Additionally, the northern burial appears to be at the same ALS elevation as the very top of the shell-rich sand layer. However, this layer does not appear level with the burial at TLS elevations.

Discussion
We encounter significant noise in a GPR dataset from a cemetery. Railroad ties that surround two graves produce ringing that overprints a anomaly of interest. SVD ranks structures in a GPR image in order of energy contribution. Because the ringing is such high amplitude and so spatially limited, it is captured in the first handful of eigenimages and, thus, easily omitted in data reconstruction.
SVD successfully removes the direct wave and the ringing to reveal overprinted diffractions. This improves the ability of the migration scheme (in this case, 2D Kirchhoff migration) to collapse single diffractions to point reflectors rather than incorporating signal from multiple features and producing erroneously large anomalies. Migrating the SVD-filtered GPR images yields more distinct and spatially constrained point reflectors. Furthermore, because signal ringing and diffractions produced by reflections are due to differing wave behavior, the ground velocity for a diffraction would only have made ringing more of a nuisance. Removing it allows for cleaner migration results.
SVD as a filtering technique is user-driven and data-dependent. One can rely on the first eigenimage to remove the direct wave, as other authors have shown. However, subsequent choices meant to remove unwanted signal will depend on (1) whether the amplitude of the reflected signal from a target of interest and "noise" are similar or (2) whether the unwanted signal exhibits the same geometry as the target signal. In archaeological settings, SVD may be successful in isolating targets such as graves and refuse pits from layered features, like shell midden or compacted ground. Again, one would need to produce the eigenimages and decide what can be filtered and what cannot.
At such a small scale (<150 square meter area), subtle relief can influence small geophysical targets' depth and shape. We show that the marked burials are not associated with a depression, as measured by ALS. Instead, they are situated beneath a subtle modern mound bounded by the railroad ties as measured by TLS. ALS data, collected in 2007, have approximately 20 cm vertical accuracy. The horizontal accuracy is 1 m. We did not know of any physical changes that would have altered the site's relief. We attribute the elevation differences to source data resolution.
The targets of interest are at approximately 20 cm higher elevation if we perform a topographic correction using TLS versus ALS. The difference is slightly larger (30 cm) for the northern confirmed burial. Under some circumstances, this is an insignificant variation. However, in archaeological applications, a notable change is often recorded at the sub-centimeter (certainly sub-meter) scale. Thus, the highest vertical resolution possible is preferred.
When presenting 3D GPR data as depth slices (no elevation applied), the elevation is inconsequential because the depth to a target is dependent on estimated ground velocities. Data interpretations from elevation slices, on the other hand, requires more careful attention in terms of the elevation source and resolution. In the example presented here, the confirmed burials have a similar vertical relationship, regardless of the elevation source. However, they are closer in depth, according to TLS versus ALS-derived elevations. Additionally, other features, like tree roots and the top of a shell-rich sand layer, are situated at different elevations. If these features were the intended target, the elevation source would influence in their depth and extent as well.

Conclusions
Geophysical interpretations are dependent on data quality and confidence in the location of anomalies relative to each other. This is especially true of GPR data from an archaeological setting where anomalies are small and shallow, and depth slices are preferred modes of presentation. Through the joint application of a subspace projection method and static topographic correction using surface elevation tomography, we show increased spatial resolution of two burials from a cemetery with two original headstones and a third feature that may or may not be associated with the burials. Terrain variations between the two elevation sources presented are small, but can be significant in the context of archaeological investigations.  Acknowledgments: Thank you to Alexandra Kulenguski, Brett Parbus, Danielle Waite, and Danielle Young for their assistance collecting field data. Thank you to Jorge González García and Kyutae (Simon) Ahn for their assistance in preparing TLS data.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of this study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

ALS
Airborne laser scanning CCAFS Cape Canaveral Air Force Station DEM Digital elevation model GPR Ground-penetrating radar GPS Global positioning system SVD Singular value decomposition TLS Terrestrial laser scanning

Appendix A. Singular Value Decomposition
Singular value decomposition (SVD) is a factorization of an m × n matrix X into orthogonal matrices in the form: where U is the matrix of column eigenvectors of XX T and V is the matrix of row eigenvectors of X T X. Σ is a diagonal matrix of singular values of X is decreasing order. A powerful feature of SVD is the ability to take high-dimensional, highly variable data, and reduce it to a lower-dimensional space that reveals substructures. We treat a GPR image as an m x n matrix where m is the number of traces, and n is the number of samples per trace. It can be written as: where r is the rank of X (the number of linearly independent traces), σ is ith singluar value of X, u i is the ith eigenvector of U, and v i is the ith eigenvector of V.
Eigenimages are linearly independent components of the GPR image. Singular values are the energy contribution factor of a given eigenimage. For example, the first (and largest) singular value corresponds to the eigenimage that contributes the most to the whole composition, which, in GPR, is he direct wave [10,15,38].