Continuous Coastal Monitoring with an Automated Terrestrial Lidar Scanner

Abstract

This paper details the collection, geo-referencing, and data processing algorithms for a fully-automated, permanently deployed terrestrial lidar system for coastal monitoring. The lidar is fixed on a 4-m structure located on a shore-backing dune in Duck, North Carolina. Each hour, the lidar collects a three-dimensional framescan of the nearshore region along with a 30-min two-dimensional linescan time series oriented directly offshore, with a linescan repetition rate of approximately 7 Hz. The data are geo-referenced each hour using a rigorous co-registration process that fits 11 fixed planes to a baseline scan to account for small platform movements, and the residual errors from the fit are used to assess the accuracy of the rectification. This process decreased the mean error (defined as the magnitude of the offset in three planes) over a two-year period by 24.41 cm relative to using a fixed rectification matrix. The automated data processing algorithm then filters and grids the data to generate a dry-beach digital elevation model (DEM) from the framescan along with hourly wave runup, hydrodynamic, and morphologic statistics from the linescan time series. The lidar has collected data semi-continuously since January 2015 (with gaps occurring while the lidar was malfunctioning or being serviced), resulting in an hourly data set spanning four years as of January 2019. Examples of data products and potential applications spanning a range of spatial and temporal scales relevant to coastal processes are discussed.

1. Introduction

Collecting data in the nearshore environment is a notoriously difficult task. Breaking waves and strong nearshore currents make nearshore sampling challenging and dangerous, and the relevant hydrodynamic and morphodynamic processes vary over a broad range of spatial and temporal scales. Traditional in situ measurement and survey techniques generally provide spatially sparse datasets, which often fail to capture the spatial variability in wave and current fields, bathymetry, and beach morphology necessary to understand and properly characterize many of these processes [1]. These data collection methods are also primarily deployed over relatively short time periods (days to months, in the case of fixed instruments) (e.g., [2]) or at a low temporal resolution (monthly to yearly, in the case of surveying) (e.g., [3]) because of financial or logistical constraints. These approaches limit the utility of these datasets to assess longer-term change (months to decades) at the time scales needed (seconds to days) to improve our physics-based understanding of coastline evolution.

To address this need, continuous remote sensing methods have been developed that can provide high resolution coastal data over large areas for long durations [4]. For example, shore-based optical observations of the surf-zone such as Argus [1,5] have been used to extract information on nearshore morphology [6,7,8]; to estimate bathymetry [9], beach slope [10], alongshore currents [11], and wave runup [12]; and to quantify shoreline change over long timescales [13]. Although the value of continuous remote sensing is clear, single-location optical systems such as Argus cannot provide highly accurate measurements of water surface and beach elevations, which are important parameters in both coastal science and management.

In contrast, airborne, mobile, or terrestrial lidar systems are capable of collecting elevation data in a wide range of coastal environments, and their use has become increasingly common in recent years [14,15,16,17]. In coastal areas, lidar scanners have been used to produce both large-scale coastal maps [18,19] as well as small-scale high resolution scans of nearshore morphology [20,21,22]. Additionally, lidar scanners are capable of providing returns off both the ground and the water surface, allowing for the assessment of beach and dune morphology (e.g., [20,21,22,23,24,25]) as well as wave transformation and water surface elevation in the swash and inner surf-zones [26,27,28,29,30].

Although the use of lidar in coastal science and monitoring has rapidly increased in the past decade, most lidar systems are deployed over a fixed time period (days to weeks) (e.g., [31,32]) or with low frequency (monthly to annually) (e.g., [20,25]) due to financial, technological, or logistical constraints. However, recent advances in scanner technology and computing capabilities, reductions in the cost of equipment and data storage, and improvements in data processing methods have extended the duration that lidar systems can feasibly be deployed and facilitated the development of automated lidar systems for geophysical research [33,34]. This paper describes the automated data collection and processing algorithms for a fixed lidar system deployed at the US Army Engineer Research and Development Center’s Field Research Facility, the data products generated by the system, and potential applications of this dataset in coastal science and monitoring efforts. This method represents a step forward in coastal observing capabilities, as beach elevation change data can now be autonomously collected at the same frequency (seconds to hours) as the driving physical forcing parameters (waves and currents).

2. Methods: Data Collection and Processing

2.1. Field Site

The US Army Engineer Research and Development Center’s Field Research Facility (FRF), part of the Coastal and Hydraulic Laboratory (CHL), is a coastal observatory with a 560 m research pier located near the town of Duck, NC, on the Atlantic coast of the United States. The facility was established in 1977 to support field investigations in the coastal sciences and has since been the site of a large number of coastal field studies (e.g., [2,35,36,37]. The FRF maintains a number of permanentlydeployed in situ instruments and remote sensors that provide continuous hydrodynamic and morphodynamic data to the coastal community. Data are publicly available and can be accessed on the CHL THREDDS data server (CHL TDS) (https://chlthredds.erdc.dren.mil/thredds/catalog/frf/catalog.html). The location of the FRF property, the FRF research pier, and the town of Duck are shown in Figure 1.

2.2. Scanner Setup and Scan Collection

Scans are collected using a Riegl terrestrial VZ-1000 lidar scanner (1550 nm laser with a 0.3 mrad beam width) housed in a Riegl VZ-1000-PH protective case. The scanner is deployed above the foredune crest at the north end of the FRF property (located at 36.1858° N, 75.7524° W and x = 47.86 m, y = 945.33 m in FRF coordinates, shown in Figure 1). The scanner is mounted at an elevation of 13.8 m in the North American Vertical Datum of 1988 (NAVD88) on the end of a mast angled 30° above the horizontal (Figure 1, inset), which is generally located between 20–40 m from the waterline. The landward end of the mast is mounted to a vertical support structure that is anchored into a concrete pad, and the seaward end rests in a support cradle which is designed to limit the movement of the scanner. The mounting point is hinged, allowing the mast to be lowered to access and service the lidar scanner. The lidar system is shore-powered, and commands to the lidar system and data from the lidar system are routed using a fiber-Ethernet path running from the scanner to the FRF facility.

The lidar scanner collects a 30-min linescan time series and a three-dimensional spatial framescan spanning 237° around the scanner each hour. In this set-up, cross-shore linescans are collected at 7.1 Hz, a pre-determined function of the peak laser pulse repetition rate (70 kHz), vertical angular resolution (0.025°), and vertical collection range (100°). The scans provide a 30-min time series of the elevation of the dry beach, swash, and inner surf-zone waves along this transect each hour. The 30-min collection period was chosen to adequately sample most wave-driven hydrodynamics within an appropriate window of assumed stationarity while still allowing time for the subsequent framescans to be collected.

Framescans are collected over a 16-min period following the linescan time series collection each hour. To ensure a relatively consistent horizontal data density along the beach, three separate framescans are collected and merged, allowing for variable azimuthal and vertical angular resolution in different areas around the scanner (Figure 2). The portions of the beach located farther from the scanner (included in Scan 1 and Scan 3, shown in Figure 2) are sampled at a higher vertical angular resolution (0.01°) than Scan 2 (0.04°). The azimuthal angular resolution is 0.05° in Scans 1 and 3 and 0.04° in Scan 2. All scans are collected with a peak laser pulse repetition rate of 70 kHz. Point densities on the beach range from 10² to 10⁴ points/m² within 100 m of the scanner and decrease with range to 10¹ points/m² far (e.g., 400 m) from the scanner (Figure 2).

The lidar scanner has been collecting scans semi-continuously since January 2015, with gaps occurring while the scanner was malfunctioning or being serviced. As of January 2019, the scanner has collected 24,573 hourly scans out of 34,680 total hours, generating a dataset that includes data from 15 storm events (defined as events with a significant wave height > 4 m in 26 m water depth).

2.3. Processing Methods

2.3.1. Geo-Referencing and Scan Alignment

All linescans and framescans are collected in the Scanner’s Own Coordinate System (SOCS) and are then transformed into a geo-referenced coordinate system used locally at the FRF. The FRF coordinate system has the positive x-axis pointing offshore (aligned with the FRF research pier, which can be seen in Figure 1 and Figure 2) and the positive y-axis pointing 18° west of true north, with NAVD88 as the vertical datum. After the hourly framescans are collected, the automated geo-referencing algorithm transforms the scans into FRF coordinates and then co-registers the scan to a baseline scan to account for any movement of the scanner between scans. The transformation matrix generated in the co-registration process is used to rectify both the linescan and framescan for that hour, and the accuracy of the co-registration is assessed. A flow chart outlining the steps in the scan geo-referencing and alignment process is shown in Figure 3 and a detailed explanation of the processing steps is provided in the paragraphs to follow.

All framescans are co-registered to a single baseline framescan (collected on 27 July 2017 using the same scanner parameters as the hourly framescans) to ensure accurate rectification and continuity in the position of fixed objects through time. Immediately following the baseline scan, fine scans of 19 cylindrical retroreflector targets (with a height of 10 cm and a diameter of 10 cm) were also collected. Five of these retroreflector targets are permanently deployed and are used in the error assessment described in Section 2.3.2. The remaining targets were temporarily deployed during the baseline scan collection. A survey of each reflector was conducted using real-time kinematic (RTK) global positioning system (GPS) surveying equipment and the measured GPS locations were converted to the local FRF coordinate system. RiScan Pro, Riegl’s operating and processing software, was used to identify the reflectors in the baseline framescan. A rigid transformation matrix RT_B was then determined using a least-squares adjustment between the reflector locations in SOCS and the measured FRF coordinates (with a standard deviation of 0.0265 m between the measured reflector centers and the reflector center coordinates in the scan). RT_B was then used to transform the baseline scan from SOCS to FRF coordinates using the equation

$$\left[\begin{array}{c}{x}_{FRF}\\ y_{FRF}\\ z_{FRF}\\ 1\end{array}\right]=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_x\\ r_{21}& r_{22}& r_{23}& t_y\\ r_{31}& r_{32}& r_{33}& t_z\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_{SOCS}\\ y_{SOCS}\\ z_{SOCS}\\ 1\end{array}\right]={\mathbf{RT}}_{\mathbf{B}}\left[\begin{array}{c}{x}_{SOCS}\\ y_{SOCS}\\ z_{SOCS}\\ 1\end{array}\right]$$

(1)

where r_ij are the rotational components and t_i are the translational components of RT_B. Once the baseline framescan was transformed to FRF coordinates, the scan was trimmed using pre-determined x-, y-, and z- bounds to include only points on 11 planar surfaces throughout the scan (black circles show location of planes, Figure 2). The planes are all parts of solid structures (including houses, buildings on the FRF property, and beams on the FRF pier) that are not expected to move. The planes were selected to be as well-distributed around the scanner as possible to ensure a well-constrained fit during least-squares minimization. These 11 planes from the baseline scan are used as the control plane surfaces in the co-registration process, described in detail in the following paragraphs.

After each subsequent hourly scan is collected, the framescan point cloud is roughly transformed to the FRF coordinate system using an initial rectification matrix RT_I before the scan is co-registered. A new RT_I is generated any time the scanner is removed from the mount for maintenance (occurring on average once every two months) to account for any minor changes in the resting location of the scanner once it is remounted. Similar to the baseline rectification method, a least-squares fit of the reflector centroids in the scan to the measured reflector locations is used to solve for RT_I. RT_I is then used to initially rectify all hourly scans collected until the scanner is again removed for maintenance.

Once an hourly scan has been initially rectified, the scan is trimmed using the same pre-determined x-, y-, and z- bounds to include only points on the 11 fixed planes. The points in the trimmed scan are conditioned to lie on the control plane surfaces (the same planes from the baseline scan, defined using the centroid and normal vector of each baseline plane) through the adjustment of the six parameters of a rigid body transform. Covariance matrices generated from the propagation of raw measurement error to the coordinates of each point in the trimmed hourly scan (following Hartzell et al. [38]) are used to weight the points in the least-squares adjustment using the equations

$$C_{xyz}=\left[\begin{array}{ccc}{\sigma }_x{}^2& {\sigma }_{xy}& {\sigma }_{xz}\\ {\sigma }_{yx}& {\sigma }_y{}^2& {\sigma }_{yz}\\ {\sigma }_{zx}& {\sigma }_{zy}& {\sigma }_z{}^2\end{array}\right]=A{C}_{ll}{{\rm A}}^{\mathrm{T}}$$

(2)

where σ are the covariances between dimensions x, y, and z; C_xyz is the unregistered point covariance matrix; C_ll is the observation covariance matrix, determined from the laser beam width and uncertainties in range and angle; and A defines a linear relationship between the observations (range, vertical angle, and azimuthal angle) and unregistered point coordinates (point coordinates in FRF coordinates prior to co-registration). Because the relationship between the observations and the unregistered point coordinates is non-linear, partial derivatives are used to estimate a linear relationship, with A defined as

$$A=\left[\begin{array}{ccc}\cos \theta \cos \varphi & -r\cos \theta \sin \varphi & -r\sin \theta \cos \varphi \\ \cos \theta \sin \varphi & r\cos \theta \cos \varphi & -r\sin \theta \sin \varphi \\ \sin \theta & 0& r\cos \theta \end{array}\right]$$

(3)

where θ is the vertical angle and ϕ is the azimuthal angle.

Using the three rotation parameters (α₁, α₂, α₃) determined from the least-squares fit of the 11 planes, a co-registration rotation matrix, R_c, is generated using the following equation:

$$R_c=\left[\begin{array}{ccc}\cos {\alpha }_2\cos {\alpha }_3& -\cos {\alpha }_2\sin {\alpha }_3& \sin {\alpha }_2\\ \cos {\alpha }_1\sin {\alpha }_3+\sin {\alpha }_1\sin {\alpha }_2\cos {\alpha }_3& \cos {\alpha }_1\cos {\alpha }_3-\sin {\alpha }_1\sin {\alpha }_2\sin {\alpha }_3& -\sin {\alpha }_1\cos {\alpha }_2\\ \sin {\alpha }_1\sin {\alpha }_3-\cos {\alpha }_1\sin {\alpha }_2\cos {\alpha }_3& \sin {\alpha }_1\cos {\alpha }_3+\cos {\alpha }_1\sin {\alpha }_2\sin {\alpha }_3& \cos {\alpha }_1\cos {\alpha }_2\end{array}\right]$$

(4)

R_c and the translation vector T_c are then combined with the initial transformation matrix RT_I to form a final transformation matrix that includes both the initial rectification and co-registration matrices, which may be expressed as

$$\left[R{T}_F\right]=\left[R_c{R}_I|T_c+T_I\right]$$

(5)

where R_c and R_I are the rotation matrices from the co-registration and initial rectification, respectively, and T_c and T_I are the translation vectors from the co-registration and initial rectification. The final transformation matrix RT_F is used to transform the hourly framescan and linescan from the SOCS coordinate system to the FRF coordinate system using eq 1. In some cases, the algorithm does not locate any points on one or multiple planes. This may be a result of a shadowing effect from some object between the scanner and the plane or low point densities due to rain or fog. In these cases, the algorithm attempts to co-register the scan with the remaining located planes, as long as at least five planes are found.

The co-registration process described above removes much of the error in scan geo-rectification induced by slight movements of the scanner between scans (quantified in detail in the following section). For example, in Figure 4a, a large positive bias relative to the baseline scan is visible in the elevation differences between an example framescan (26 July 2017 17:00 UTC) rectified with RT_I on the terrain to the north of and behind the scanner (red) and a large negative bias is visible to the south of the scanner (blue) in areas which should be relatively stable. In contrast, once the scan is rectified with RT_F (Figure 4b), the offsets are significantly reduced (with the assessment plane errors, discussed in Section 2.3.2, decreasing from 0.44 m in the scan rectified with RT_I to 0.0041 m in the scan rectified with RT_F). In Figure 4b, the visible colors are now primarily a result of actual beach change and returns off the water surface, which are removed during the filtering process (discussed in Section 2.3.3).

2.3.2. Co-registration Error Assessment

Several assessment metrics were developed to evaluate the accuracy of the co-registration methodology presented above and to identify scans with high geo-rectification errors after co-registration, which generally occur during poor weather conditions (e.g., rain, wind, and fog). This automated process utilizes several empirically-determined and site-specific error thresholds which have been determined based on an analysis of 3.5 years of data and a comparison between error metrics. Although these numbers are site- and use-specific, similar thresholds could be developed for other locations.

First, the co-registration error is assessed using the residual standard errors in translation Tσ, calculated using the normal matrix N and references variance matrix S₀² determined during the least-squares fit of the 11 co-registration planes. The reference variance matrix S₀² is found using

$$S_0{}^2=\frac{K\text{'}WK}{f}$$

(6)

where K is the vector of residuals; W is the weighting matrix determined from the unregistered point covariance matrix C_xyz and the x-, y-, and z-coordinates of all points located on the co-registration planes; and f is the degrees of freedom. The residual standard errors e_residual are calculated for each hourly scan using

$$e_{residual}=\sqrt{S_0{}^2*diag(N^{-1})}$$

(7)

If the magnitude of the residual standard error in translation (defined as X = $\sqrt{T{\sigma }_x^2+T{\sigma }_y^2+T{\sigma }_z^2}$) is over 0.4 cm, the scan is flagged as part of the quality control process. This threshold only includes the residual standard error in translation and does not incorporate the residual errors in rotation, and therefore is not representative of the actual error in the scan. However, the residual errors in rotation have been found to follow similar trends, and this threshold has been found to be a reliable indicator of a poor least-squares fit and thus high geo-rectification error (Figure 5b,c).

Second, the accuracy of the co-registration is assessed using three additional planes in the scan which are independent from the 11 planes used in the least-squares fit. Each of the three planes is aligned in FRF coordinates such that they each are two-dimensional with a fixed coordinate in the third FRF dimension. The assessment planes include a support beam on the FRF pier (a y-z plane with a fixed x-coordinate), the front of a structure on the FRF property (an x-z plane with a fixed y-coordinate), and the pavement surface in a cul-de-sac in the neighborhood to the north of the scanner (an x-y plane with an approximately fixed z-coordinate). The scan offset in each direction is determined using these planes, which provides a more physically intuitive indicator of co-registration error. As an example, the offset in the y-direction on the x-z error assessment plane between the baseline scan and an example framescan rectified using RT_I (before co-registration) and RT_F (after co-registration) can be seen in Figure 4c,d. If the magnitude of the offset in the three planes exceeds 10 cm, the scan is flagged.

Third, the accuracy of the scan co-registration is assessed using the location of five permanent reflectors in the scan, ranging in location from 23 to 108 m from the scanner (red circles, Figure 2). The location of the centroid of each reflector face is found by trimming the scan down to 1 m³ around the surveyed location of each reflector. The points with the highest reflectance (3.5% of the total points in each cube) are considered reflector points, and the centroid of these points is determined. The centroid of each reflector face is then compared to the centroid location of the same reflector in the baseline scan. If the root-mean-square (RMS) error of the five reflector centroids exceeds 10 cm, the scan is flagged.

During the period from 15 November 2015 through 15 November 2017 (two years), a total of 11,434 out of 12,242 scans (93.4% of all scans) were co-registered using the plane-matching method described above with error metrics within the acceptable range. Prior to co-registration, the magnitude of the assessment plane offsets in the scans that were eventually successfully co-registered ranged from 0.0244 cm to 132 cm, with a mean of 26.2 cm and a standard deviation of 18.0 cm. After co-registration, the magnitude of the assessment plane offsets ranged from 0.0248 cm to 9.86 cm in the same scans, with a mean of 1.77 cm and a standard deviation of 1.38 cm. A subset of the co-registration errors (for the period from 1 December 2016 through 7 January 2017) can be seen in Figure 5, including the magnitude of the assessment plane offsets (Figure 5a), the residual standard errors in translation (Figure 5b), and the residual standard errors in rotation (Figure 5c). A diurnal spike in error values was observed in the rectified scans prior to co-registration (red line, Figure 5a), potentially due to the daily thermal expansion of the base on which the lidar is mounted. The co-registration process was found to remove most of this daily error increase (compare blue and red line, Figure 5a).

2.3.3. Scan Filtering

Once the scans have been rectified in space, all unwanted points are removed before the scan is interpolated onto a regular grid. The linescans are conducted at a high temporal resolution that allows the surf-zone water surface (particularly shoaling and breaking waves) to be resolved. However, the slow azimuthal rotation rate during the framescan collection process results in little usable water or wave information. For this reason, all returns from the water surface in the framescan are removed. To do this, the cross-shore position of the instantaneous water line is identified at each alongshore location and all data offshore of this location are removed. Water tends to have significantly lower reflectance values than dry beach (Figure 6a), creating a relatively continuous and coherent change in reflectance at the water-beach line. The waterline is identified in the scan by applying a Canny edge detection algorithm using the reflectance values of each point. The algorithm is first applied using a high reflectance threshold and then run iteratively with decreasing thresholds until a continuous waterline is identified along the beach. During this iterative process, the waterline identification is guided by constraining the distance offshore and onshore that the waterline can move between adjacent cells. The location of the waterline identified by the filtering algorithm can be seen in the example filtered digital elevation model (DEM) shown in Figure 6c. This line matches the location of the rapid change in reflectance at the waterline seen in Figure 6a.

Once water is removed from the framescan, a second filtering process is applied to the scan to remove noise, objects, and dune vegetation. First, all points with reflectance lower than −25 dB are removed from the scan, which removes most returns from particles in the air. Next, a cloth simulation filter [39] is used to separate ground and non-ground points. The cloth simulation filter inverts the point cloud and fits a flexible “cloth” with a defined rigidness and a 50 cm window to the underside of the point cloud (see Zhang et al. [39] for filter details). This surface is then considered the ground, and all points over 10 cm off this surface are removed. This filter removes vegetation; people, objects, and animals on the beach; and any remaining noise (blue dots, Figure 6c,d). The points remaining after the filtering process (red circles, Figure 6c,d) are considered ground points and are used to generate the bare earth DEM.

A separate filtering process is applied to the linescan data to remove noise, objects, and sea spray, as well as multiple-reflection points. The linescans are first separated into a dry beach region and a swash/inner surf-zone region, and the two regions are filtered separately. The edge of the dry beach is determined by finding the cross-shore location where the variance of reflectance values increases markedly (because points are alternatingly returning off wet beach and water, which have different reflectance characteristics), and this point is considered the start of the swash zone. Within the dry beach region, the filter removes all points over 10 cm from the elevation of the maximum point density at each cross-shore location, which removes vegetation and temporary objects from the beach.

The swash zone and water surface in the linescans are filtered using a series of empirically defined conditions tuned for this location. First, all points are removed that appear at elevations below the ground or water surface at each cross-shore location. These points are primarily the result of an emitted pulse reflecting off multiple surfaces (often off the smooth water surface and then off the face of an incoming wave), and thus taking longer to reflect back to the scanner. In the swash zone, these points are removed by finding the elevation of the wet beach and removing points below this. Offshore of the swash zone, these points are removed by assuming that any points with an elevation more than a specified distance below all other data at the same cross-shore location over the 30-min collection period are multiple-reflection points. Next, sea spray and noise above the water are removed using the assumption that each consecutive return should be returning from a cross-shore location farther from the scanner than the previous return. Only points that have both a cross-shore location closer to the scanner than the previous point as well as a difference in elevation of at least 10 cm from the previous point are removed to avoid removing the curling edge of breaking waves. Finally, the location where the point density falls below 100 (over a 10 cm interval) during the 30-min collection is found, and all points seaward of this are removed. Examples of filtered and unfiltered linescan time series can be seen in Figure 7.

2.3.4. Point Cloud Gridding

The rectified and filtered framescan point clouds are next interpolated onto a regular grid using a nearest neighbor binning technique. The regular grid spans 500 m in the alongshore direction and 130 m in the cross-shore direction with 1 m spatial resolution. This grid is produced operationally every hour; however, finer-scale grids can be produced if needed at any time.

The linescans are interpolated onto a 10 cm cross-shore grid extending from x = 42.9 m to x = 197.8 m in FRF coordinates using a linear interpolation algorithm. The interpolation algorithm does not interpolate over a gap spanning a cross-shore distance greater than 10 m or an elevation jump of greater than 5 m. A discussion of point density and potential shadowing/interpolation issues is included in Brodie et al. [28] and is therefore omitted here. The final product is a continuous sea-surface from each individual linescan from which hydrodynamic statistics can be calculated.

3. Results: Data Products

The gridded, filtered framescans and linescan time series are then automatically processed to quantify several hydrodynamic and morphodynamic parameters, as outlined below. These parameters are publicly available in NetCDF files stored on the CHL TDS (https://chlthredds.erdc.dren.mil/thredds/catalog/frf/catalog.html).

3.1. Runup

The lidar linescan time series provides information about both the hourly changes in beach morphology (Section 3.3) as well as the extent of the wave runup. Wave runup is a function of both the time-averaged water level at the shoreline as well as the fluctuations about that mean (i.e., swash) [40] and is a primary driver of beach and dune erosion [41]. Wave runup can be quantified by the lidar [26] by generating a 30-min time series of the position (cross-shore and alongshore coordinates) and elevation (NAVD88) of the instantaneous waterline at approximately 7 Hz. While differences in reflectance are useful for determining the approximate beach-water interface for filtering (see Section 2.3.3), the similarity in reflectance between wet-sand and water (not shown) prevents the use of reflectance for automatically determining the instantaneous runup at the higher sampling rate of the linescans. Instead, changes in elevation through time at each cross-shore location are used to separate water and beach data following Brodie et al. [28] to generate the runup time-series. Points seaward of this location are classified as water and points landward of this location are classified as land. This automated process can contain errors in the location of the runup line. For this reason, the runup time-series is manually checked to assess the quality and to correct its location if needed. Once manually checked, a quality control flag is changed on the CHL TDS to indicate that data is appropriate for scientific use by the community. Improving the automated runup finder and removing this manual step is a future research goal. In addition to the runup time-series, the horizontal and vertical R2% is determined for each 30-min linescan. R2% is defined as the cross-shore location and elevation exceeded 2% of the time by wave runup events and provides a statistical representation of the maximum onshore extent of wave runup during the time series [12].

3.2. Hydrodynamics

The hourly hydrodynamic measurements determined from the scans include water surface elevation statistics, spectral analyses, and third-order moments of the elevation time-series. Only cross-shore locations with high returns (defined as locations with returns in at least 75% of the time series) are included in the analysis. These analyses are conducted from the runup line (discussed above) every 10 cm to the offshore extent of high returns, meaning both the onshore and offshore extent of hydrodynamic information varies based on wave and tide conditions. The water surface statistics include the mean water level and the total significant wave height H_s (where $H_s=4\sigma$ and σ is the standard deviation), both of which are determined at each cross-shore location.

A spectral analysis of the water surface elevation time series at each cross-shore location is also conducted. The spectral analysis provides the power spectral density estimate of the wave energy, calculated using Welch’s method with 4.8 min windows, 75% overlap between windows, and band-averaging over three frequency bands, resulting in 36 degrees of freedom. From the wave spectrum, the significant wave height H_s and mean period T_m are determined in both the infragravity band (f < 0.04 Hz) and in the wind sea and swell band (0.04 Hz < f < 0.5 Hz). H_s at each band is defined as $H_s=4\sqrt{m_0}$ where m₀ is the zeroth moment (determined by integrating the power spectral density over the given frequency band). T_m is determined by finding the frequency associated with the energy centroid in each band (with T_m = 1/f_m).

Finally, the wave skewness and asymmetry are calculated at each cross-shore location. Wave skewness refers to the vertical asymmetry in the waveform (i.e., differences between the shape of the crest and trough), and is calculated using the equation

$$s=\frac{\mathrm{E}{(\eta -\mu )}^3}{{\sigma }^3}$$

(8)

where η is the water surface elevation, µ is the mean of η, σ is the standard deviation of η, and E(η − µ)³ represents the expected value of (η − µ)³. The wave asymmetry, or horizontal asymmetry in the waveform (i.e., differences between the face and back of a wave), is defined as the skewness of the imaginary component of the Hilbert transform of the surface elevation time series [42].

While the hourly hydrodynamic statistics and spectral analyses described above are provided across the total extent of the data, due to space constraints, the full time series of water surface elevation are only provided at a subset of cross-shore locations on CHL TDS. These virtual wave gauges are spaced every 10 m from x = 80 m to x = 140 m in FRF coordinates. Quality control metrics indicate the percentage of the total time series with no measurements as well as the median duration of the data gaps.

3.3. Morphology

The framescan processing results in an hourly bare earth DEM (Section 2.3.4). The co-registration error metrics indicating the quality of the scan co-registration are reported along with each DEM. The hourly morphological parameters determined from the linescans include the mean and standard deviation of the foreshore elevation and the foreshore slope. The foreshore elevation mean/standard deviation is defined as the mean/standard deviation of the elevation of the beach over the 30-min collection at each cross-shore location, in 10 cm increments from the scanner. The foreshore slope in the swash zone is defined as the slope of the best-fit linear trend to the beach elevation profile in the region defined by the mean horizontal swash position plus or minus two standard deviations in the swash [12].

4. Discussion

4.1. Applications

Advances in remote sensing technologies have expanded the types, durations, and resolutions of datasets that can be collected in coastal zones and thus the types of questions that can be addressed using these data. The method outlined in this paper results in a dataset that can be used in hydrodynamic or morphodynamic analyses on scales of seconds to years and includes continuous data during both storms and quiescent periods. Examples of potential applications of the continuous lidar framescan and linescan datasets are outlined below.

4.1.1. Long (Seasonal to Annual) Time Scales

Potential analyses of low frequency processes (i.e., tidal cycles to years) using this dataset include assessments of seasonal erosion and accretion trends, beach cusp formation, and the evolution of foreshore slope, among others. Furthermore, the high spatial resolution and large spatial extent of the present dataset allow for the assessment of the spatial variability of these processes across and along the beach. An example of the cross-shore variability in morphodynamic processes can be seen in the elevation time series shown in Figure 8b, taken from a fixed alongshore location (y = 950 m) and from cross-shore locations ranging from x = 55 to 95 m in FRF coordinates, with blue colors located higher on the beach and red colors located lower on the beach. The figure shows the relative stability of the upper beach time series, with change primarily occurring during large wave events (with H_s in 8 m water depth shown in Figure 8a). The lower beach is much more dynamic, with change occurring over much a wider range of temporal and spatial scales.

4.1.2. Event Scales

In addition to analyses of processes occurring over long time scales, the lidar dataset described in this paper could also be used to conduct analyses on the time scale of single storms. Storms and other extreme events often result in rapid changes in beach and dune morphology. Because of the difficulty collecting data during large storms, few observations of surf-zone hydrodynamic or morphodynamic processes during extreme storms exist [27,31]. For this reason, models used to simulate coastal change due to storms often use parameterizations calibrated using only pre- and post-storm bathymetry and topography [43,44,45,46]. The framescan and linescan datasets described in this paper include data collected during several large storm events, including Hurricanes Joaquin, Matthew, Jose, and Maria (Figure 9), among others. Additionally, because of the continuous nature of the dataset, logistics concerns associated with rapid deployment are removed—the instrument is already in place when the storm arrives. Analyses of beach and dune evolution during storms could be used to help validate or improve numerical models simulating extreme storm processes and impacts.

For example, coastal engineers and scientists often strive to predict total water level during storms, which is a sum of tide, surge, and wave-driven R2% [12]. Measured wave-driven R2% (red dots, Figure 9a) was calculated during three storms with offshore significant wave heights exceeding 4.5 m: Hurricanes Joaquin (left panel), Jose (middle panel), and Maria (right panel). Interestingly, while offshore wave heights were similar between the three storms, the wave-driven R2% was significantly higher during Joaquin (with a maximum of 3.75 m during Joaquin relative to 2.52 m and 2.32 m during Jose and Maria, respectively). This resulted in higher levels of erosion at the base of the dune during Joaquin, visible in both the hourly beach elevations transects (left panel, Figure 9b) and in the cross-shore elevation change timestacks (left panel, Figure 9c). The transects and elevation change timestacks (Figure 9b,c) show a spatial variability in erosion and accretion patterns during the storm, with accretion occurring on the upper beach face during Jose and Maria (middle and right panels, Figure 9b,c). The runup statistics determined from the linescan dataset could be used to determine the validity and applicability of commonly used wave runup parameterizations (e.g., [12]) in a range of morphologic and hydrodynamic conditions (e.g., [47]).

4.1.3. Short (Wave-by-Wave) Time Scales

The high temporal and spatial resolutions of the linescan dataset provide additional opportunities for analyses on time scales down to individual waves. As waves propagate toward the coast and begin to shoal, their profiles evolve from generally sinusoidal in deep water to positive skewed (i.e., sharper crests and flatter troughs) in intermediate water depths, and eventually to positive asymmetric (i.e., steeper front faces and more gently sloping rear faces) near and in the surf-zone [42]. These wave profile changes result in asymmetric orbital velocities and accelerations, which can affect both the magnitude and direction of cross-shore sediment transport [48]. Past work has shown that pressure- and lidar-based measurements of wave skewness and asymmetry are qualitatively consistent with one another, suggesting that lidar measurements can be used to analyze wave profile evolution across the surf-zone [28]. The linescan dataset described in this paper provides wave profile information at a high spatial and temporal resolution across the inner surf-zone, which could be used to further investigate the evolution of shoaling and breaking wave shapes near the shoreline [29,30]. Figure 10 shows an example of the evolving profile of two waves (separated by 20 min in time) as they propagate through the surf-zone. In the top panel, the evolution of the waveform from positive skewed farther offshore to positive asymmetric (right to left) as the wave approaches the shoreline is clearly visible. In this case, the wave evolves from positive skewed near 125 m in the cross-shore to positive asymmetric only 25 m onshore in roughly 4 s, highlighting the rapid evolution on timescales less than a wave period and spatial scales shorter than a wavelength. In the bottom panel, lidar data show two separate bores merging as they approach the shoreline (blue to red from x = 95 m to x = 70 m, Figure 10b). Bore-bore capture has been shown to play a role in controlling wave runup patterns and other surf-zone processes [49]. The long duration of the dataset provides an opportunity to analyze inner-surf wave shapes in a range of wave conditions and over a range of surf-zone and beach morphologies.

In addition to the hourly analyses of wave runup statistics and morphologic change discussed in Section 4.1.2, the linescan dataset can also be used to extract information on morphologic change in the swash zone induced by individual runup events. The swash zone is the site of a significant amount of total surf-zone sediment transport [50]. Although a number of studies have been conducted that focus on sediment transport in the swash zone (e.g., [51,52,53,54,55,56,57,58,59,60]), few have quantified the morphological change in the dune toe induced by individual swash events. Figure 11a shows the swash depth through time during a 4.5 min period starting on 4 October 2015 18:17:45 UTC. Of these runup events, 17 reached an elevation above the dune toe (shown with a magenta line in Figure 11b). The sediment flux into or out of the dune was calculated for these 17 events (Figure 11c) along with the cumulative volume change (Figure 11d), which demonstrate the variability of the sediment flux induced by individual events during this period. Additionally, the cross-shore transects of beach elevation change shown in Figure 11b show the spatial variability of erosion and accretion patterns across the beach face, with an area of accretion (from x = 55 m to x = 63 m) occurring between the generally erosive areas at the dune toe (x < 55 m) and on the lower beach face (x > 63 m).

4.2. Error Propagation

In order to fully interpret and exploit the lidar data collected using the method outlined above, steps must be taken to understand and quantify the scale of errors across the scan. The co-registration method described above provides the additional benefit of generating transformation parameter covariance information during the least-squares fit, which can be propagated along with estimated errors in the lidar sensor range, angle, and laser beam width into the lidar point cloud [38] to provide spatially variable error estimates for the final data. To evaluate whether the observed offsets on our assessment planes could be explained by these uncertainties, we propagated the error through the scan following the General Law of Propagation of Variance (GLOPOV) outlined in eq 2 to calculate a standard error in each assessment plane centroid in the direction of interest at the 95% confidence interval. Because the relationship between the observations and registered coordinates is non-linear, the matrix A is defined using a Taylor Series linear approximation of the relationship between these parameters. In this formulation, the spatially variable error estimates are strongly dependent on both the range from the lidar as well as the orientation of the surface relative to the lidar location. For an example period from 20 November 2016, 00:00 UTC to 10 December 2016, 00:00 UTC, the calculated standard errors in the assessment plane centroid coordinate components were found to bound the observed plane offsets at the 95% confidence interval in 243 out of the 259 available scans in the x-direction, 242 out of 259 scans in the y-direction, and 65 out of 259 scans in the z-direction. Although many of the observed scan offsets exceeded the calculated standard error in the assessment plane centroid in the z-direction, the mean exceedance was less than 2 mm.

This observed error exceedance may occur because the error propagation does not compensate for errors in lidar measurements due to atmospheric effects or slight motions of the scanner during the scan. Additionally, observed errors may be slightly higher than the propagated error due to the potentially incorrect assumptions that the assessment plane surface is flat (zero surface roughness), completely stationary, and perfectly aligned in FRF coordinates. The ability to calculate spatially variable uncertainty estimates that scale with actual errors is useful for downstream uses of the data products such as assimilation into numerical models or including error bounds in volume change estimates. Future work will continue to investigate these results to determine when and why the propagated error bounds or does not bound the observed error, and how optimistic or conservative the estimate is.

5. Conclusions

This paper describes the data collection and processing algorithms for a permanently-deployed coastal lidar system located at the US Army Corps of Engineers Field Research Facility (FRF) in Duck, NC. The automated processing algorithm includes a rigorous co-registration and error assessment process that ensures that the resulting data are properly geo-rectified. The unprecedented resolution and duration of this dataset resulting from the method described in this paper provide opportunities for analyses of coastal morphodynamic processes at a wide range of temporal and spatial scales. Examples of potential analyses are given at time scales ranging from seconds to years. On the time scales of individual waves, potential analyses of the linescan dataset include assessments of wave asymmetry and skewness over a range of conditions and surf-zone bathymetries. On longer time scales, the framescan dataset could be used to analyze seasonal erosion and accretion trends as well as beach cusp and scarp formation, and to monitor beach volume changes following a nourishment. Finally, both framescans and linescans could be used to assess hydrodynamic and morphodynamic processes during extreme events such as storms, including spatially variable accretion and erosion patterns across the beach face as well as dune volume changes per swash events.