Submeter Resolution Surface Rupture Topography From Legacy Aerial Photographs—A Test Case From the 1992 Landers Earthquake

The 1992 Mw 7.3 Landers earthquake in the Mojave Desert (California) provided exceptional observations of surface faulting in a large, continental earthquake. The U. S. Geological Survey obtained nadir angle, overlapping aerial photographs at 1:6,000 scale for the entire ∼ 85 km rupture length. Recent advances in Structure from Motion photogrammetry allow for archival photographic data sets such as these to be reprocessed, generating digital topography that can be reanalyzed quantitatively in a way that was not previously possible. In this proof‐of‐concept study, we generated a georectified, ∼ 10 points/m 2 topographic point cloud over nearly the entire Landers rupture length and a higher‐resolution ∼ 40 points/m 2 point cloud over a smaller ( ∼ 5 km) rupture segment along the Emerson fault. We estimated the accuracy and explore the utility of our point cloud in two tests. First, we observe close geometric agreement (average closest point distance 2.1 cm and standard deviation 14.0 cm) between our point cloud and a 2008 terrestrial lidar survey of the Galway Lake Road site on the Emerson fault. Second, we made 173 vertical offset measurements within a small, structurally complex pull‐apart basin, also on the Emerson fault, and find visual and statistical similarity with 21 local field measurements. These two tests demonstrate that point clouds generated from legacy aerial surveys and georeferenced using free Google Earth and National Elevation Dataset imagery are geometrically accurate and can be used to densify geomorphic offset measurements even along well‐studied surface ruptures. Applied to other historical events, such measurements could provide new insights into earthquake rupture processes.

From south to north, the Landers earthquake sequentially ruptured the Johnson Valley, Kickapoo, Homestead Valley, Emerson, and Camp Rock faults (Figure 2a). Few of these faults were well characterized before the earthquake, and the Kickapoo fault had not been recognized as active. The intricate and well-exposed surface ruptures were surveyed aerially by the United States Geological Survey (USGS) 2 days after the earthquake, principally to assist with field-based mapping (Arrowsmith & Rhodes, 1994;Aydin & Du, 1995;Kaneda & Rockwell, 2009;McGill & Rubin, 1999;Sieh et al., 1993;Spotila & Sieh, 1995 (Sieh et al., 1993). These multimeter dextral offsets were often accompanied by vertical offsets of ∼1 m or more. Dextral slip on the northernmost rupture segment, the Camp Rock fault, was mostly limited to <1 m.
Starting around a decade ago, airborne lidar topography has transformed our means of detecting subtle tectonic landforms, extracting fault offset measurements, determining fault surface geometries, and geomorphic mapping of slip rate and trenching sites (e.g., Meigs, 2013;Zielke et al., 2015). For such applications, lidar has been used extensively across Southern California (e.g., Blisniuk et al., 2012;DeLong et al., 2010;Frankel et al., 2007;Hudnut et al., 2002;Oskin et al., 2007;Oskin et al., 2008;Salisbury et al., 2012;Zielke et al., 2010). Most major known faults in California now have airborne lidar coverage (e.g., Bevis et al., 2005;Crosby et al., 2013), but conspicuously and perhaps surprisingly, the majority of the Landers rupture has never been scanned ( Figure 1). The absence of high-resolution digital imagery has limited further analysis of its surface rupture characteristics and may explain the near absence of Late Quaternary slip rate estimates along its constitutive faults (Oskin et al., 2008).
Recent advances in optical image analysis provide a new means of generating high-resolution topography along the Landers rupture using the existing photographic imagery. Structure from motion (SfM) is a method that uses two-dimensional optical images with varying viewing geometries to create three-dimensional point clouds (irregular sets of external surface coordinates) that, in sparsely vegetated regions, approximate the topography (James & Robson, 2012;Snavely et al., 2008;Westoby et al., 2012). Unlike classical stereophotogrammetry, image locations and camera parameters need not be specified beforehand but instead are determined along with point cloud coordinates in a single "bundle adjustment" calculation. This technique enables reconstruction of landscapes at high spatial resolution using archival photosets lacking detailed metadata (Bakker & Lane, 2017;Derrien et al., 2015;Gomez et al., 2015;Midgley & Tonkin, 2017).
This study presents a proof-of-concept approach for generating useful topographic data along historical earthquake ruptures from archival aerial photographs, using the Landers earthquake as a test case. We construct a ∼10 points/m 2 point cloud and submeter resolution digital elevation model (DEM) from the overlapping USGS aerial images. The only archived geographical information for these photographs are approximate corner and center points of the photographs, which provide limited and imprecise information for properly registering them and their derived products. However, we show that modern, freely available geographical information provides an alternative means of georeferencing archival photosets. Next, we use published field survey data to evaluate our point clouds, demonstrating that high-resolution digital topography is geometrically accurate and can potentially yield new fault offset measurements along even well-studied earthquake ruptures. Similar overlapping aerial photosets likely exist for numerous modern 10.1029/2019EA000651 and historic earthquakes, and we show that these data can be used to create high-resolution topographic models that rival airborne lidar in resolution and accuracy.

Data
We utilized a USGS aerial photographic survey flown over almost the entire ∼85 km Landers rupture trace in the late morning to early afternoon of June 30, just 2 days after the earthquake. Black-and-white, approximately nadir angle photographs were collected at a scale of 1:6,000. Our analysis is based on 365 high-resolution digital scans of the nine-by-nine inch negatives of the 1:6,000 scale images. Though corner and center points of each photograph are archived, their accuracy cannot be verified; because this information is not essential to the SfM processing, we disregarded it.
We draw upon two other published data sets in this study. First, we used a subset of a small (∼200 m by ∼250 m, or ∼0.04 km 2 ) terrestrial lidar scanning (TLS) survey collected in 2008 at Galway Lake Road on the central Emerson fault (https://doi.org/10.5069/G91Z4298; Figure 2c). These data are used to determine whether our point clouds are geometrically robust. The TLS point cloud was scanned from 11 ground positions using a tripod-mounted Riegl LPM 321 terrestrial laser scanner. A total of 8.8·10 6 points are reported with point densities of up to ∼ 38 · 10 4 points/m 2 and averaging ∼231 points/m 2 . Eighteen GPS ground control points (GCPs) allowed the survey to be reprojected into the 2008 Universal Transverse Mercator (UTM) Zone 11 N NAD83 coordinate system with a <10 cm uncertainty in network adjustments. Vertical measurements are referenced to the WGS84 ellipsoid.
Second, we draw upon coseismic offsets mapped on the ground by McGill and Rubin (1999) along the central Emerson fault. We used these field data to assess whether or not our topographic point clouds can yield robust vertical offset measurements. Aided by the 1:6,000 USGS aerial photography, McGill and Rubin (1999) used a steel tape to make 60 offset measurements across a ∼5.6 km section of the main fault zone and another 200 across secondary structures. McGill and Rubin (1999) focused on fractures with ≥10 cm of slip, though they measured smaller offsets that occurred on structures with larger slips in other localities.
An equivalent, pre-1992 topographic model of the rupture area might allow us to compute the three-dimensional surface deformation field for the Landers earthquake and more accurately assess fault slip and coseismic surface displacements (e.g., Lajoie et al., 2019;Nissen et al., 2014). Toward this end, we attempted to produce a preevent topographic point cloud using 51 black and white images collected in nine aerial surveys between August 1947 and June 1975 at scales ranging from 1:30,000 to 1:136,000 (USGS "Earth Explorer" data portal, http://earthexplorer.usgs.gov/). Unfortunately, the noise level of the resulting point cloud was greater than the expected multimeter deformation signal and was therefore useless for any quantitative purposes. However, we note that pre-earthquake and postearthquake photographic imagery were successfully used to determine the horizontal coseismic displacement field (Milliner et al., 2015).

Point Cloud Generation
We used Structure-from-Motion (SfM) software to reconstruct the postearthquake surface topography from the scanned aerial photographs, generating high-resolution ground surface point clouds with coregistered texture (color). SfM computes the structure of a scene by triangulating the positions of features recognized in multiple images viewed from different angles and distances. It makes use of a feature recognition algorithm, the Scale Invariant Feature Transform (Lowe, 2004;Snavely et al., 2008), that can accommodate irregular or unknown image acquisition geometries as well as changing illumination and partial occlusion. Importantly, "structure" refers to both the geometry of the imaged feature (topography, for aerial surveys in sparsely vegetated landscapes) and camera parameters including location relative to the scene, orientation, and camera lens focal length. Consequently, SfM can be used even when no metadata exist, making it ideally suited for use with legacy data sets where such information is crude or lacking.
Following the work of Johnson et al. (2014), we constructed SfM point clouds using the commercial software Agisoft Photoscan Pro, version 1.3.0 (for a list of alternative software, including some that are freely available, see Bemis et al., 2014). While the workflow within Photoscan is largely automated, the lack of geographic image metadata requires that ground control be added manually. We approximated 75 GCPs throughout the coverage area of the 1:6,000 aerial survey by coupling Google Earth imagery with the National Elevation Dataset (NED) to determine x, , and z coordinates of recognizable features-mostly low-lying bushes and small boulders-located in broad, flat areas. Ideally, all geographic information would be hosted in a single source, but NED image resolution is insufficient to locate small features despite having much higher-resolution elevation data than Google Earth. Google Earth incorporates smoothed elevations from the 30 m Shuttle Radar Topography Mission; the NED elevations are resolved at 1/3 arc second (about 10 m at the latitude of California). Conveniently, geographic coordinates for both data sets are referenced to the WGS-84 ellipsoid. We determined UTM easting and northing coordinates of GCPs using Google Earth, and then obtained corresponding elevation values using a bulk point query to the NED. The purpose of choosing GCPs located in broad, flat areas is to minimize uncertainties introduced by elevation data smoothing, but we also ensured that collectively the GCPs encompass a wide range of elevations. We estimated lateral uncertainties for the GCPs of ∼1 m based on our confidence in visually identifying targets across data sets; the vertical uncertainty is more difficult to assess, and likely exceeds ∼1 m.
Ground control can be imported into Photoscan Pro without visual reference; GPS measurements, for example, can be uploaded and used to manually rescale, rotate, and translate a point cloud. While we acknowledge that this might provide more accurate ground control in areas where such acquisitions are possible, we elected to use only georeferencing data that are readily available for any geographic locality, thus maximizing the applicability of this proof-of-concept test.
Owing to large file sizes (averaging ∼60 MB), the 365 individual 1:6,000 images were divided into eight subsets that were merged into a single point cloud after final processing. Each photo subset was processed with the feature matching algorithm set at Photoscan's "medium accuracy" setting, because the higher accuracy settings increased computation time and became prohibitive for such large data sets. Accuracy in this case refers to the number, not the quality, of matched features. Similarly, when building the dense point clouds, we elected the "high quality" rather than "ultrahigh quality" processing setting, again reducing the resolution of the resulting point cloud along with the computation time, but not the quality of modeled points. The raw point cloud was denoised using the Photoscan "gradual selection" tool on each of the sparse point clouds. This tool calculates the reprojection error between the estimated and real location of a point in a pixel and is a metric of the quality of the camera alignment and subsequent feature extraction. The cleaned, full-rupture length point cloud has an average density of ∼10.3 points/m 2 (Figure 2b).
To test higher-resolution capabilities of this method, we generated a smaller point cloud with a subset of 17 images that encompasses the southern study area (the orange polygon in Figures 2b and 2c). For this region, we placed 5 GCPs with ∼1 m accuracy and processed the data using Photoscan's feature matching algorithm set to "high accuracy" and dense cloud quality set to "ultrahigh". The resulting point cloud has ∼329 million points and covers ∼8 km 2 with an average point density of ∼40.1 points/m 2 , a fourfold increase in spatial resolution over the full-rupture point cloud (Figure 4).

Point Cloud Accuracy
A scientifically useful topographic data set must properly represent the geometry of the original ground surface, with both shape and scale accurately rendered. Ideally, geographic referencing is also accurate, though meaningful measurements can be made without absolute geographic location provided that the internal geometry is correct. In this section, we assessed the geometric accuracy of the full-rupture point cloud by comparing a small subset to a colocated 2008 TLS survey. This comparison serves as a proxy for the accuracy of the point cloud as a whole, though we acknowledge that variations in topography and image coverage will affect local point cloud quality. The TLS point uncertainties are <10 cm, so for the purpose of this assessment, we treated the TLS scan (with 18 GPS GCPs) as a "ground truth" representation of the target surface. Accuracy was quantified by computing closest point distances from the SfM point cloud to the TLS point cloud; very small residual distances between topographically complex surfaces require that the SfM point cloud is both properly scaled and geometrically accurate. For perfectly planar surfaces, however, very small distances between point clouds can be achieved even with large horizontal registration errors and scale discrepancies, such that the distance between registered point clouds is no longer relevant metric of similarity. Fortunately, most of the test area contains short wavelength, decimeter-to meter-scale relief, and so this is not a concern of ours.

Methods
We calculated distances between the SfM and TLS point clouds using the open-source point cloud modeling software CloudCompare (https://www.danielgm.net/cc/). We trimmed the SfM point cloud to match the area of the TLS survey ( Figure 3a). The TLS coverage is sparse in the very steepest areas, and so we additionally removed points from the SfM cloud in these areas, leaving ∼ 7.3 · 10 5 points in the SfM cloud. To remove any effects of geographic misregistration between the two surveys, we first finely registered the two point clouds to each other using an implementation of the iterative closest point (ICP) algorithm in CloudCompare (the "fine registration" tool). ICP uses iterative rigid body transformations (rotations and translations of the SfM point cloud to minimize distances to closest points within the TLS cloud-in this case requiring 2.3, 4.1, and −10.5 m x, , and z translations, respectively. All rotations were negligible, in no case exceeding ∼ 3 · 10 −3 degrees. After registration, we computed distances from our SfM points to a meshed surface of the TLS survey using the CloudCompare "cloud-to-mesh" tool.

Results
Following ICP closest point alignment, we compare the two surveys using a map of cloud-to-mesh distance from the registered SfM point cloud to the meshed TLS survey (Figure 3b). The average cloud-to-mesh distance is 2.1 cm with a standard deviation of ∼14.0 cm. These centimeter level values demonstrate that the multimeter translations needed to align the clouds with ICP mostly represent differences in geographic projections for ground control between the two surveys, as opposed to internal geometric inaccuracies in the SfM point cloud. Positive cloud-to-mesh values highlight the scarp and stream channels, likely areas of erosion and redistribution of sediment in the 16 years between surveys. Indeed, repeat topographic surveys and ground stereoscopic photography at this locality in the 11 months following the earthquake found as much as 50 cm of lateral scarp knick point migration in the first year alone, due largely to heavy rainfall in the winter of 1992-1993 (Arrowsmith & Rhodes, 1994). We can therefore attribute some of the morphological difference between the 1992 SfM cloud and the 2008 TLS cloud to erosion, rather than to geometrical errors in either point cloud. The minimal cloud-to-mesh distances and ability to observe erosional patterns confirm that our SfM point cloud is a properly scaled and geometrically accurate representation of the true postearthquake surface topography.

Vertical Offset Measurements
We also assessed the utility of our point clouds by attempting to replicate field measurements made after the earthquake. For this purpose, we used a smaller, higher-resolution point cloud, and chose a ∼700 m-long stretch of the central Emerson fault rich in structural complexity and with dense field measurements as a test area (Figures 2c, 4, and 5). At our selected location the northward propagating right lateral rupture stepped right to form a small (up to ∼100 m wide) pull-apart basin. We used 21 field offset measurements from ∼200 m south to ∼200 m north of the pull-apart basin, which show right lateral coseismic offsets across the main rupture zone of between 9 and 434 cm and a maximum vertical separation of 100 cm, eastside up (McGill & Rubin, 1999). While field data for the entire rupture have been digitized and compiled by Liu et al. (2003), we found substantial geographic misregistration between their measurements and our surfaces, so instead of  Rubin, 1999). A two-sided Kolmogorov-Smirnov test cannot reject the null hypothesis that these two data sets are from the same continuous distribution at the 5% significance level, indicating that point cloud-derived measurements are statistically similar to field measurements.
using their data we manually digitized field measurement locations from McGill and Rubin (1999). However, McGill and Rubin (1999) noted high spatial variability in their slip measurements, and so even small errors in location might cause significant differences in measured vertical offset. Consequently, we measured offset at 10 m intervals along the fault and statistically tested whether the two sets of measurements are from the same population.

Methods
We measured vertical fault offsets from 68 evenly spaced swath profiles through our SfM point cloud along the pull-apart basin (Figure 5a). We approximated the fault using three generalized segments, with the central segment cutting diagonally across the pull-apart basin. Profiles were oriented perpendicular to the closest fault segment and were centered at 10 m increments along the generalized fault trace. Each swath profile is 150 m long and 1 m wide. Offsets were measured by regressing linear least squares fits through point cloud profiles on either side of scarps, extrapolating the best fit lines to the vertical center point of the scarp, and computing the vertical separation between intersection points. We developed our own MATLAB code to semiautomate this process, allowing us to dictate scarp locations and the lateral intervals over which lines are fit to the point cloud, and giving us the discretion to select only bare ground.
Previous studies have found a correlation between vertical offset in the 1992 earthquake and tectonically uplifted alluvial terraces suggestive of repeat vertical offsets (e.g., Bull, 1996;McGill & Rubin, 1999). The eastern side of the pull-apart basin is flanked by such terraces, that have retreated away from the fault. Since we seek to compare SfM-derived scarp heights with those measured in the field after the 1992 earthquake, we must ensure that our SfM profiles capture single-event (1992) scarps and avoid older surfaces. However, scarp degradation studies at this site suggest that scarps as large as those generated in the 1992 earthquake would take about 10 kyr to erode to an unidentifiable state (Arrowsmith & Rhodes, 1994), and nearby paleoseismic trenching dated the most recent faulting event to ∼9 ka (Rubin & Sieh, 1997). Therefore, the alluvial surface within and immediately surrounding the pull-apart basin was likely close to planar immediately before the earthquake. Accordingly, we consider that our profiles display only vertical displacements associated with the 1992 Landers event. Figure 5a shows the locations, magnitudes, and polarities of both field and point cloud vertical offset measurements superimposed on a hillshaded SfM topographic basemap; Figure 5b shows these data projected onto a N-S profile. To the 21 field measurements we added 173 point cloud-derived offsets that capture deformation patterns both within and adjacent to the pull-apart basin (supporting information, Table S1). Point cloud offset measurements clearly delineate the intricate network of normal faults mapped by field teams and visible in the hillshaded SfM DEM (Figure 4b). Profiles advancing from north to south across the pull-apart basin highlight the large oblique-normal displacements along the western margin of the basin that were not captured by field measurements (Figure 6). Some of the profiles (e.g., C-C ′ ) also show apparent block rotations within the basin, with tilting of the basin floor toward the closest bounding faults.

Results
The field and point cloud-derived scarp offsets are in good agreement, given uncertainties in the colocation of the two measurement sets (Figure 5d). A histogram of vertical offsets for both data sets shows a bimodal distribution of point cloud offsets, with peaks around −20 cm (i.e., westside up) and +20 cm (eastside up) ( Figure 5c). Figures 5b and 5c hint at a tendency toward slightly greater field offsets than point cloud offsets. This discrepancy might indicate that field measurements were made where the scarps were clearest, likely coinciding with local peaks in surface slip. A two-sided Kolmogorov-Smirnov test is unable to reject the null hypothesis that the two data sets are from the same continuous distribution (or identical distributions) at the 5% significance level, confirming in a statistical sense that our point cloud profiling method yields offsets comparable to field measurements.

10.1029/2019EA000651
Our surface rupture offset mapping also confirms a pattern in the field measurements-that the scarp immediately south of the pull-apart basin is westside up. This behavior only persists for ∼300 m south of the pull-apart according to field measurements (McGill & Rubin, 1999), reverting back to eastside up farther south, near the edge of our test area. Our results confirm that vertical offsets are significantly larger north of the pull-apart basin (averaging ∼50 cm) than south of it (averaging ∼10 cm), increasing across the step-over as the rupture propagated northward.

Estimation of Coseismic Deformation
With the first order the accuracy of our topographic point cloud confirmed, we can use it to explore the structure and deformation of the pull-apart basin. In section 5.1, we explained that the subplanar active alluvial surface in the area of the pull-apart basin was probably underformed before the 1992 earthquake (Arrowsmith & Rhodes, 1994;Bull, 1996;McGill & Rubin, 1999;Rubin & Sieh, 1997). Disregarding low-relief stream channels, this surface can be approximated as a gently dipping (∼ 2 • ) plane. This approximation allows us to calculate the deflection of the SfM point cloud from that idealized plane, and thus estimate the coseismic vertical deformation.
We first masked out the older, raised alluvial surfaces, as well as a thin, highly deformed strip ∼10 m either side of the principal 1992 fault trace. The remaining point cloud defines a mostly planar surface that we fitted with a plane dipping 2 • toward an azimuth of 038 • . The point cloud deviation from this surface was computed for the complete (unmasked) point cloud using the cloud-to-mesh calculation tool in CloudCompare ( Figure 7). The magnitude of this deviation along the fault is low to 0 at the southern edge of the test area, increasing northward with westside-up sense of offset at the southern edge of the pull-apart basin, and the fault emerges from the northern end of the basin with consistent eastside-up offset of up to 135 cm. Away from the fault trace, surface deflections reflect the nonplanar nature of the surface, the planar approximation being an oversimplification that is only useful for qualitative assessments. Nonetheless, the large vertical deflections that define the pull-apart basin are robust and show that the basin is deepest slightly westward from a bisecting line through the long axis of the basin, closer to the large westside-up offsets than the smaller eastside-up offsets along the basin flanks in Figure 5. We note that this pattern is contrary to the one implied by the original field data set, in which measured fault offsets are much greater along the eastern wall of the pull-apart.

Conclusions
In this paper, we show that it is possible to make accurate and scientifically relevant topographic point clouds by applying SfM to sets of legacy air photos georeferenced with freely available remotely sensed data. We use the 1992 Landers rupture-one of few major faults in California currently lacking in airborne lidar coverage-to illustrate. At the Galway Lake Road site on the central Emerson fault, differences between our SfM point cloud and a GPS-georeferenced terrestrial lidar survey are mostly in the order of a few centimeters, reflecting a combination of survey errors and geomorphic change between the two surveys. At the nearby pull-apart basin site, we are able to reproduce field measurements of vertical offset within statistical significance and also densify these data substantially over what has been published previously. We demonstrate that point clouds derived from historic air photos can yield new offset measurements that potentially add to our understanding of even the best studied historic earthquakes.