INTRODUCTION

Hurricane Sandy, known unofficially as “Superstorm Sandy,” made landfall as an intense post-tropical cyclone on the U.S. East Coast near Brigantine, New Jersey, on October 29, 2012 (Halverson and Rabenhorst, 2013; NOAA, 2013). Factors contributing to the devastating impact of the storm included its very large diameter, its impact angle, and the fact that its landfall in the New Jersey–New York region coincided with large astronomical tides to produce exacerbated storm tides (Forbes et al., 2014; Hall and Sobel, 2013). Immediate impacts of the storm included at least 147 deaths, $50 billion in damages, and extensive coastal erosion in New Jersey, New York, and other mid-Atlantic states (Blake et al., 2013; NOAA, 2013). Long-term ecological impacts of the storm are still being assessed. Coinciding with the location of the center of the cyclone at the time of landfall on the U.S. East Coast, the Barnegat Bay estuary was heavily impacted by Sandy. The bay experienced ~2 m of storm surge and extensive damage, dune erosion, massive property damage, and deposition of marine debris in the estuary (Blake et al., 2013; Miselis et al., 2013). Due to the extent of damage, Barnegat Bay has become a focal point for a number of studies related to Hurricane Sandy. Coastal zone management offices are interested in assessing the effects of Hurricane Sandy on benthic habitats, particularly seagrass habitats, in Barnegat Bay.

Seagrasses are important for the health of estuarine systems, as they provide habitat for fish and shellfish species, reduce sediment erosion and currents, and deliver nutrients from the estuary (Orth and van Montfrans, 1987; Zimmermann, 2003). In Barnegat Bay, greater nutrient loading has led to estuarine-wide declines in seagrass populations and greater abundance of nuisance and non-native macroalgal species (Fertig, Kennish, and Sakowicz, 2013; Hunchak-Kariouk and Nicholson; 2001; Kennish, Fertig, and Lathrop, 2012; Kennish, Fertig, and Sakowicz; 2011; Seitzinger, Styles, and Pilling, 2001; Weiben and Baker, 2009). It has been suggested that seagrass beds are expanding into shallow subtidal areas due to reduced light availability by phytoplankton, macroalgal blooms, and epiphytic overgrowth (Lathrop and Haag, 2011; Lathrop et al., 2001). These areas are vulnerable to damage caused by Hurricane Sandy, as the storm surges resulted in ocean water washing over parts of the barrier islands, creating strong currents and sediment deposition that could exacerbate seagrass decline by smothering or fragmenting existing seagrass beds. Lidar-derived data products such as those used in this study can be useful to coastal managers for rapid assessment of critical habitats following storm events as they require little additional postprocessing beyond standard lidar and imagery products.

Traditionally, seagrass beds have been mapped using RGB aerial imagery and field surveys (e.g., Lathrop et al., 2014; Lathrop, Montesano, and Haag, 2006; Macleod and Congalton, 1998). Aerial and satellite imagery acquired in clear water conditions enable the operator to resolve patches of vegetation in shallow waters (Moore et al., 2000). However, these sensors are passive systems that depend on solar illumination conditions, as well as water clarity (Pe'eri et al., 2016). Field surveys are beneficial as they provide on-the-ground inspection that enables detailed assessment of vegetation and species. While these surveys are extremely useful for reference data acquisition, they only provide information at discrete, sampled locations, making it difficult to assess habitat change over large spatial extents.

Bathymetric or topobathymetric lidar provides an alternative survey tool to assess temporal and/or spatial changes in habitat as a result of disturbance (e.g., storms or invasive species). Lidar is an active survey technology that is independent of the ambient illumination. A major advancement in bathymetric lidar over the past decade has been the development of radiometric calibration algorithms and procedures that enable generation of seafloor reflectance images (Macon et al., 2008; Tuell et al., 2005; Wang and Philpot, 2007). These lidar-derived seafloor reflectance images have proven useful in coastal and benthic habitat mapping (Chust et al., 2013; Costa, Battista, and Pittman, 2009), and can be produced for some areas too deep or turbid to map with passive multi- or hyperspectral imagery acquired from an aircraft or satellite.

To date, research on seafloor reflectance mapping with lidar has focused primarily on conventional bathymetric lidar systems, which use relatively high transmit pulse power, wide beam divergence, wide receiver field of view, and low pulse rate to achieve bathymetric measurements in waters up to 2–3 Secchi depths (Chust et al., 2013; Tuell et al., 2005). However, there is an emerging class of topobathymetric lidar systems that are closer in design characteristics to conventional topographic lidar than to conventional bathymetric systems. These systems use only a single, green laser with low transmit pulse energy, narrow beam divergence, narrow receiver field of view (FOV), and high measurement rates. These topobathymetric lidar systems do not offer the same depth penetration as the conventional bathymetric lidar systems, but they enable cost-effective, high-density bathymetry to be acquired in shallow waters. This type of general design was pioneered with the Experimental Advanced Airborne Research Lidar (EAARL) system (Brock et al., 2004; Heidemann et al., 2012) and further advanced by the U.S. Geological Survey (USGS) EAARL-B team (Wright, 2014).

NOAA's post-Sandy mapping efforts included data acquisition with a new, commercial topobathymetric lidar system, the Riegl VQ-820-G. This dataset did not include recorded waveforms (i.e., digitized samples of the entire backscattered signal for each transmitted pulse). The lack of recorded waveforms and the differences in system design between the VQ-820-G and conventional, large-FOV bathymetric lidar systems make it impossible to directly apply reflectance mapping procedures that have been used in other studies. However, two waveform features, termed reflectance and pulse deviation, were recorded in real time and made available as point attributes in output point clouds. Leveraging these features, the objectives of this study were to

The motivation for avoiding rigorous radiometric calibration stems from the fact that 1) for hurricane response applications, it is highly beneficial to have short turnaround times for generating habitat change maps, at least for preliminary analysis; and 2) coastal zone management offices, which constitute one of the primary intended user groups for this work, typically lack the signal processing expertise and software to perform extensive radiometric pre-processing.

Riegl Lidar Waveform Features

In Riegl's “V-line” laser scanners (all of which, other than the VQ-820-G and newer VQ-880-G, are designed for topographic or terrestrial applications), waveform processing is done in real time, and full waveforms are not always available to end users for further processing. However, as noted above, two waveform features are computed during acquisition: reflectance, ρ, and pulse deviation, δ (Pfennigbauer and Ullrich, 2010). In the lidar community, the term “intensity” is used to refer to the stored value of the peak amplitude for a given return and, in the case of a linear detector, it is proportional to the received optical power. Reflectance, as defined by Riegl, is arguably better termed “twice-normalized intensity.” The first normalization is performed by taking the ratio of the peak return amplitude (or echo power), P_echo, to the minimum-detectable signal level, P_min, with the output expressed in decibels:

The second normalization is designed to correct for the strong range-dependence of A_dB. Specifically, the reflectance is normalized by that of a white reference target, oriented perpendicular to the direction of propagation of the laser beam, at the same range, R:

The values of P_min and A_dB(R) are system specific and are computed and tabulated by Riegl in their calibration laboratory (Riegl, 2014).

The second waveform feature, pulse deviation, δ, measures the difference in area under the curve (numerical integration) between the received pulse, y[n], and a stored reference pulse, p[n]:

In Equation 3, the variable n denotes the discrete time index or sample number. As with P_min and A_dB(R), the reference pulse or “system response function,” p[n], is system specific and measured by Riegl at their factory. Both reflectance and pulse deviation can be stored as point attributes in the American Society of Photogrammetry and Remote Sensing (ASPRS) LAS file format for point cloud data (Graham, 2005; Samberg, 2007) using LAS v1.4 ExtraBytes variable length records (VLRs) (Riegl, 2014).

While the motivation for avoiding rigorous radiometric calibration in this study was stated above, the theoretical basis for doing so can be explained as follows. From the bathymetric lidar equation given in various forms in a number of studies employing data from conventional, large-FOV bathymetric lidar systems (e.g., Kopilevich et al., 2005; Tuell et al., 2005), received signal power varies linearly with seafloor reflectance, whereas there is exponential attenuation with depth. Hence, an area that contains brighter or darker values than the surrounding areas in Riegl relative reflectance datasets could well be shallower or deeper, rather than of a different benthic habitat type (e.g., eelgrass, as opposed to sand). However, from analysis of previous studies, existing habitat maps, and our own field data, areas of eelgrass habitat in Barnegat Bay correspond to a narrow range of shallow depths. Therefore, it may be possible to use the reflectance data without any depth correction in eelgrass habitat mapping. The underlying assumption is that, as a first-order approximation, the waveform features in shallow-water areas vary primarily as a function of habitat type. Furthermore, preliminary visual inspection of gridded maps of waveform feature data revealed that texture and shape are more useful in segmenting different benthic habitats than absolute brightness. For example, eelgrass and macroalgae may have similar brightness values in a lidar reflectance image (even if rigorous radiometric calibration has been performed), but they were observed in visual analysis to have quite different texture in the relative reflectance images. Finally, even though a depth correction is not explicitly performed, the bathymetry layer is available for use in the rule set development within this OBIA procedure.

Study Site

The study site in this work is the Barnegat Inlet flood tidal delta complex (Figure 1). During the storm, this area was extensively flooded with concurrent strong currents, making its effects on benthic habitats of interest to investigate. The site is composed of multiple shoals and is bifurcated by several channels that appear to migrate frequently. Tidal flow through the inlet is strongly affected by shoaling in and around the flood tidal delta (Kennish, 2001). In order to stabilize the inlet, jetties were originally constructed in 1939–1940 and have been modified or reconstructed several times. In addition, regular dredging is required to maintain a navigable waterway. As a result of the jetties and other anthropogenic alterations in the area, inlet hydraulics and sediment loading, particularly at the inlet opening, has been strongly effected (Seabergh, Cialone, and McCormick, 2003). The delta is primarily composed of medium sand with coarse shell debris and some gravel lining the channels (Psusty and Silveira, 2009), which can facilitate habitat expansion across the inlet. The dominant seagrass species in this subsite is eelgrass (Zostera marina). The tidal range of 0.95 m pushes high salinity water through the inlet. Mean tidal currents at Barnegat Inlet are 1.1 m/s during flood and 1.3 m/s during ebb. Salinity in this part of the bay tends to range from 19% to 30%, with lower salinities at the mouths of rivers and creeks and greater salinities at the inlets.

Figure 1.

Study site (indicated by green rectangle) comprising the Barnegat Inlet flood tidal delta complex overlaid on Esri Online World Imagery Basemap (right), as well as the general location of the project site along the U.S. eastern seaboard (left).

A variety of submerged aquatic vegetation (SAV) species can be found in the study site, including the eelgrass (Zostera marina) (Loveland et al., 1984). The eelgrass is typically found at depths of 1 m or less on the flats around the inlets and along the backside of the barrier beach and in Manahawkin Bay (Figure 1, right). The overall composition varies greatly by season, but the consistently dominant seaweed (macroalgae) species are sea lettuce (Ulva lactuca), red seaweed (Gracilaria tikvahiae), green fleece (Codium fragile), banded weed (Ceramium fastigiatum), and red seaweed (Agardhiella subulata). Barnegat Bay estuary has been designated by the U.S. Environmental Protection Agency (EPA) as a National Estuary Program (Kennish, 2001) since 1995 and is regularly monitored by the estuary program.

Airborne Data Collection

The airborne data for this experiment were collected by NOAA's National Geodetic Survey (NGS) with a Riegl VQ-820-G in September 2013. The airborne survey platform was a NOAA DeHaviland Twin Otter. In addition to the lidar data, Applanix DSS digital aerial imagery was acquired for the project sites. Airborne data acquisition specifications are summarized in Table 1.

Table 1.

Data acquisition parameters for the September 2013 flights

Reference Data Collection

Reference data were collected by a field team in October 2013. These data included GPS coordinates and underwater camera imagery for a number of benthic habitats in and around the Barnegat Inlet flood tidal delta complex. The GPS equipment consisted of inexpensive, non-survey-grade, but ruggedized and waterproof receivers, including a Garmin eTrex 20 and WAAS-enabled EverMore SA320 USB Marine GPS receiver. Due to the use of L1-only GPS with an estimated positional accuracy of ~3 m, habitat patches of at least 10 m² were surveyed, and samples were not taken near the edge of any habitat patch. These procedures were enforced to minimize the effects of positioning uncertainty on the classification accuracy assessment. Since reference data acquisition was performed in very shallow waters (≤1.5 m), data were collected primarily from kayaks or by foot, wading through waters up to chest deep (Figure 2). Benthic habitats in the reference dataset were classified as sand, mixed sand and macroalgae, sparse eelgrass, and dense eelgrass.

Figure 2.

Field data acquisition from kayaks in Barnegat Bay.

In general, acquisition of reference data (“ground truth”) is significantly more challenging and costly for subaqueous habitat sites than terrestrial coastal sites. The primary factors in the time increase are the lower data collection rates and greater number of weather days. The main factors in the cost increase are the travel costs (due to longer field campaigns involving a greater number of personnel) and boat rental, when required. These factors can make it infeasible to acquire the number of samples per class typically recommended in the published literature on classification accuracy assessments, e.g. Congalton (1991). Nevertheless, reference data in 33–38 sites were collected for each of the following classes: eelgrass dense, mixed sand and macroalgae, and sand.

METHODS

The workflow used in conducting this study is illustrated graphically in Figure 3. The top-level inputs included the lidar bathymetry, Riegl waveform features, RGB imagery, and expert domain knowledge of the marine ecologist on the project team. The desired outputs included the final habitat map and results of the classification accuracy assessment.

Figure 3.

Full workflow used in this research. Dashed outlines indicate optional steps that were considered but not determined to be necessary in this work.

As noted earlier, one of the specific research goals of this study was to avoid performing rigorous radiometric calibration and to use the Riegl waveform features essentially “as is.” That is to say, rigorous radiometric calibration that involves correcting for attenuation in the water column and other effects was not applied. However, a small amount of preprocessing was needed to facilitate subsequent analysis by removing seamlines and large differences in contrast between adjacent swaths. This was performed through a histogram normalization, computed as follows in the case of the reflectance data:

Two potential enhancements to the pre-processing procedure were considered, but not applied, in this work. The first was to apply an incidence angle correction. Because the Riegl VQ-820-G uses a 20° forward tilt angle to maintain a nearly constant incidence angle on the water surface, reasonable results were obtained in this study without this correction, thereby reducing processing time and costs. The second potential enhancement was to apply further seamline removal in the ERDAS Imagine Fourier Transform Editor. The procedure entails computing a Fourier transform of the mosaicked waveform features (e.g., reflectance), then applying a notch filter in the Fourier domain to remove any remaining seamlines, and finally, computing the inverse transform (Parrish et al., 2014). In this work, it was unnecessary to apply this additional seamline removal in the Fourier domain, as the pre-processing steps described above did a satisfactory job of removing seamlines and other artifacts.

The next step was to produce the benthic habitat map for the study site. The input to the classification procedure consisted of 0.15-m resolution true color orthophotos from the DSS imagery, lidar-derived raster images with a 0.30-m grid size, and gridded reflectance and pulse shape deviation layers, following the preprocessing described above. The desired output consisted of a thematic map with the following five classes: 1) Low Elevation (meaning areas below the laser extinction depth), 2) Sand, 3) Dense Eelgrass, 4) Sparse Eelgrass, and 5) Mixed Macroalgae and Sand.

Classification was performed using a rule-based expert system within an OBIA software environment. The expert system was built using the Cognition Network Language (CNL), the programming language implemented in the eCognition software platform (Trimble Navigation Limited, Westminster, CO, USA). Rule-based expert systems within eCognition are referred to as “rule sets.” eCognition version 9.0.1 was used for rule set design, development, and deployment. The rule set consisted of the following main components: 1) raster processing, 2) segmentation, 3) primitive object classification, 4) primary object classification, 5) context-based refinement, 6) minimum mapping unit correction, and 7) export. Rule set development was accomplished by pairing an experienced OBIA analyst with researchers who had both domain expertise in benthic habitat and expert knowledge of the area. By collectively looking at the data, the OBIA analyst translated the knowledge of researchers into the CNL rule set following a framework developed by O'Neil-Dunne et al. (2012) for extracting features from lidar and imagery in terrestrial landscapes.

In the raster-processing step, percent slope was derived for the lidar elevation layer to create a new percent slope raster using the “surface calculation” algorithm. Three separate segmentation routines were applied in the segmentation step. In the first, those areas with an ellipsoid elevation below −40.0 m (NAD83(2011)) were segmented and classified as “Low Elevation” using the “multi-threshold” algorithm. This threshold was empirically selected as corresponding to the extinction depth of the lidar in the AOI. These “Low Elevation” areas were then excluded from all further processing. In the next segmentation routine, coarse level objects were generated. The “multiresolution segmentation” algorithm was run on all the unclassified portions of the data (i.e., not “Low Elevation”). The red, green, and blue bands for the orthophotos were used for segmentation along with the elevation, slope, and reflectance layers derived from the lidar. Each one of these layers was equally weighted. Determining the appropriate segmentation settings was an iterative process in which the settings were adjusted until the objects were as large as possible without visually appearing to cross habitat types. A second segmentation was performed to generate smaller sub-objects. The layers used in the segmentation process were the same as the larger objects, but the scale of these objects was set so as to represent individual, homogeneous features in the landscape such as a patch of algae or eelgrass.

In the third step, the small sub-objects were classified into three categories: 1) bright (generally sand), 2) dark in the imagery and lidar (primarily eelgrass with some algae), and 3) dark in the imagery alone (primarily algae with some eelgrass). The classifications were carried out using simple thresholds of imagery brightness (an average brightness of all three bands) and lidar reflectance. Using the classification information present in the sub-objects, the larger parent objects were classified based on the relative area of each of the three classes in the sub-objects layer. The parent object classification was refined in the fifth step by using the spatial relationships between objects. For example, single patches of eelgrass surrounded by mixed macroalgae and sand were reevaluated and, if warranted, reclassified. In the sixth step those objects smaller than the desired minimum mapping unit (500 m²) were dissolved into the surrounding object with the largest shared border. The classification was exported to a shapefile in the seventh and final step. The rule set developed in this work is being made publicly available via ScholarsArchive@OSU, Oregon State University's digital service for gathering, indexing, making available, and storing scholarly work.

A subset of the habitat classification map generated in eCognition is shown in Figure 4. The final classification map for the full project site, overlaid on an Esri World Imagery basemap, is shown in Figure 5.

Figure 4.

Subset of final habitat classification map for Barnegat Bay (right) and corresponding location in Google Earth imagery (left). The area shown is a ~0.7 km² site centered on Kite Island, a small, circular island, within the Barnegat Inlet flood tidal delta complex. The color scheme in this classification map is: pink = land; dark green = dense eelgrass; light green = sparse eelgrass; orange = sand; and red = water deeper than the lidar extinction depth.

Figure 5.

Final benthic habitat map for full project site.

Classification Accuracy Assessment

The habitat map was assessed by overlaying the reference data described in the Data Section as a point shapefile on the final habitat map (Figure 5) in ArcGIS and compiling an error matrix. The results of this classification accuracy assessment with all four classes are summarized in Table 2.

Table 2.

Error matrix using all four classes.

An issue with the results presented in the error matrix in Table 2 is that the “eelgrass sparse” class contained only 9 samples in the reference dataset. The reference data for this class corresponded to areas in the ground truth dataset that contained small patches of eelgrass within a larger sandflat. Hence, the 7 eelgrass sparse points that were misclassified as sand were located in an area that was, in fact, predominantly sand. The eelgrass in this location were short and occurred in small patches. This indicates that the eelgrass beds were likely expanding into this area, which is consistent with the literature that states that eelgrass beds are moving into shallower environments due to water quality (i.e., turbidity issues, as noted in Lathrop et al., 2001). Therefore, the four-class classification scheme might not be well suited for areas in which eelgrass is in the process of recolonization.

Because the small sample size for the sparse eelgrass class calls into question the statistical validity of the results, the sparse and dense eelgrass classes were collapsed into a single “eelgrass” class and the results were re-tallied (Table 3).

Table 3.

Error matrix with “eelgrass sparse” and “eelgrass dense” classes merged.

Notwithstanding the inherent challenges of acquiring reference data for subaqueous habitats, all of the classes in the revised classification scheme have >30 samples in the reference dataset. The overall accuracy for this assessment is 85%. The User's and Producer's accuracies for each class are listed in Table 4.

Table 4.

Producer's and User's accuracies for each class.

DISCUSSION AND CONCLUSIONS

The results of this classification accuracy assessment are quite encouraging. The overall classification accuracy of 85% is sufficient to support meaningful change analysis, especially when lidar surveys are available for pre- and post-storm periods. In the case of the current Sandy project, the high classification accuracy for dense eelgrass class was important for direct assessment of changes to the SAV and indirect impact to the environment (e.g., water conditions). Furthermore, the use of OBIA was found to be efficient, as the rule set development took only 7 hours, and execution of the rule set within eCognition to extract the benthic habitat classes took less than 25 minutes.

The OBIA approach to image classification used in this study precluded a direct comparison of results with vs. without inclusion of the reflectance and pulse shape deviation layers. Modifying the set of input layers would necessitate a complete redesign of the expert system, and it would be impossible to quantify the specific contribution of the data to differences in results, since the analyst plays a role. However, the OBIA analyst and ecologist who led the development of the rule set both noted, qualitatively, the importance of the reflectance layer. The value of reflectance was also evidenced by the fact that, of the 18 segmentation and classification algorithms included in the OBIA expert system, all 18 made use of lidar reflectance. The pulse deviation layer was found to be of comparatively little benefit for the benthic habitat classification in this study, contributing to only one of the algorithms.

The results indicated that there is confusion in the classification between sand (i.e., pure sand) and vegetation mixed with sand: namely, sparse eelgrass mixed macroalgae and sand. It is also difficult for an expert in the field to visually classify a mixed eelgrass/sand habitat into sparse eelgrass class or sand class. This is because the sparse eelgrass locations within the study site correspond to areas within sandflats in which eelgrass is recolonizing. Depending on the state of recolonization, it can be somewhat subjective tocategorize a site as sand or sparse eelgrass. It would be possible to partially alleviate this issue in future work using more precise class definitions (e.g., if between 30% and 60% of the habitat patch is eelgrass, then it belongs to the sparse eelgrass class). From an ecological perspective, the ability to identify any percentage of eelgrass is significant (Short, 1992).

One of the most important outputs of this work, beyond the mapping results, is the rule set developed in eCognition. To date, this rule set has been tested only on data from the Riegl VQ-820-G and a single study site. It will likely have to be modified for other sites and lidar systems. Nevertheless, this study result may serve as a useful starting point for future work involving mapping eelgrass and other seagrass habitats with topobathymetric lidar and gridded waveform features. The rule set is being made publicly-available to facilitate such studies.

The project team is currently (2015) extending this work to the USGS EAARL-B system to assess seagrass habitat change in Barnegat Bay resulting from Hurricane Sandy. Pre-Sandy EAARL-B data were acquired in the bay within three days of the storm's landfall in New Jersey. Post-Sandy acquisition started just three days after the storm and concluded within a week. Furthermore, as the EAARL-B records a complete return waveform for every pulse and receiver channel, it is possible to calculate similar waveform features to the ones investigated in this study, as well as additional features, such as pulse width and numerical integral.