1. Introduction
Knowledge of spatial snowpack properties in mountain regions, especially snow water equivalent (SWE), is critical for accurate assessment and forecasting of snowmelt timing (Reference Luce, Tarboton and CooleyLuce and others, 1998), snowmelt volume (Reference Elder, Dozier and MichaelsenElder and others, 1991) and avalanche hazard (Reference Conway and AbrahamsonConway and Abrahamson, 1984; Reference Birkeland, Hansen and BrownBirkeland and others, 1995), for initialization of synoptic and global-scale weather and climate models (Reference ListonListon, 1999; Reference Groisman, Davies, Jones, Pomeroy, Walker and HohamGroisman and Davies, 2001) and for investigations of ecologic dynamics and biogeochemical cycling (Reference Brooks and WilliamsBrooks and Williams, 1999; Reference JonesJones, 1999). Snowfall and wind interact with terrain and vegetation to create highly variable patterns of snow accumulation (e.g. Reference Winstral and MarksWinstral and Marks, 2002; Reference Grünewald, Schirmer, Mott and LehningGrünewald and others, 2010). These complex interactions produce a snow cover that is challenging to sample and model (Reference Elder, Dozier and MichaelsenElder and others, 1991). The seasonal snow system and its spatial distribution at multiple scales is coupled to hydrologic, atmospheric and biogeochemical systems through dynamic forcing of runoff characteristics, heat and energy fluxes, soil moisture distributions and growing season duration (Reference Jones, Pomeroy, Walker and HohamJones and others, 2001; Reference Brooks, Grogan, Templer, Groffman, Öquist and SchimelBrooks and others, 2011), greatly influencing energy, water and biogeochemical cycling in mountain and Earth surface systems.
Manual sampling of snow depth is expensive, time-consuming, potentially dangerous to field crews and disturbs the snowpack, influencing subsequent measurements. Avalanche starting zone depths and runout volumes are particularly difficult to sample, presenting an obvious need for remote-sensing technologies. Further, the spatial intervals over which snow depth can be feasibly measured with manual techniques are limited to spacings and extents that likely do not capture the natural variability at the slope or basin scale (e.g. Reference Elder, Dozier and MichaelsenElder and others, 1991; Reference BlöschlBlöschl, 1999; Reference Hiemstra, Liston and ReinersHiemstra and others, 2006). More recently, GPS interferometry has been useful in measuring snow depth to within 10cm (Reference Gutmann, Larson, Williams, Nievinski and ZavorotnyGutmann and others, 2012). These measurements could supplement the existing networks of snow pillows and snow courses in the western USA, but they are likewise limited in spatial extent.
Lidar (light detection and ranging) is a remote-sensing tool with the ability to retrieve precision target positions at high spatial resolutions in rough terrain and forested regions (Reference Lefsky, Cohen, Parker and HardingLefsky and others, 2002). Typically, lidar acquisitions are conducted using either airborne or ground-based instruments, and target positions are geolocated by coupling the lidar system with a high-precision GPS (ground-based) or GPS/inertial measurement unit (IMU) (airborne) system. Differencing co-registered lidar maps from two dates (snow-off from snow-on) allows the calculation of snow depth with sub-decimeter vertical uncertainty and high horizontal spatial resolutions over spatial extents compatible with basin-scale hydrologic needs and for slope-scale investigations (Fig. 1; e.g. Reference Deems, Fassnacht and ElderDeems and others, 2006; Reference Mott, Schirmer and LehningMott and others, 2011). The spatial resolution and coverage, repeatability and sub-canopy mapping capabilities of airborne and ground-based lidar offer a powerful contribution to research-oriented and operational snow hydrology and avalanche science in mountain regions.
Calculation of snow depth from lidar data requires two co-registered data collections, one each for snow-free and snow-covered dates, followed by differencing the snow surface and bare-ground elevations (Reference Miller, Elder, Cline, Davis and OchsMiller and others, 2003; Reference Hopkinson, Sitar, Chasmer and TreltzHopkinson and others, 2004; Reference Deems, Fassnacht and ElderDeems and others, 2006). Lidar-derived digital elevation models (DEMs) have accuracies as great as 10 cm root-mean-square error (RMSE) even in densely forested areas (Reference Kraus and PfeiferKraus and Pfeifer, 1998; Reference Reutebuch, McGaughey, Andersen and CarsonReutebuch and others, 2003). The snow-depth calculation procedure involves the subtraction of two surface maps plus interactions of the laser light with the snow surface. Additionally, snow-depth mapping in mountain terrain involves consideration of laser-scanning geometry relative to steep slopes and with the potential for dramatic variations in aircraft flying height and corresponding variation in ground-point spacing. These factors, if not accounted for, produce the potential to inflate horizontal and vertical uncertainties in lidar-derived snow-depth measurements compared with common accuracy estimates derived for flat snow-free surfaces.
The science of lidar mapping is evolving and the body of work concerning suitability and error sources for natural resource applications continues to grow. However, the available sensors, proprietary processing techniques and data formats are not standardized, making an understanding of each instrument, processing step and range of potential error sources critical for successful application of lidar mapping for snow science interests. This paper provides an overview of lidar measurement techniques, system parameters and error sources and magnitudes involved in snow-depth mapping using airborne and ground-based lidar, concluding with recommendations for successful employment of this powerful technology for scientific and operational snow hydrology and avalanche science needs.
2. Lidar Measurement and Snow-Depth Calculation Techniques
2.1. Lidar range measurement and system types
Lidar is an active ranging instrument and measures target distance based on the time-of-flight principle, i.e. the elapsed time between transmitted and return laser signals and energy. The position of the sensor is established by way of differential GPS (DGPS) triangulation. Additionally for airborne sensors, platform orientation (roll, pitch, yaw) is determined via an IMU (Fig. 2). Once the geometry of the platform and scanner is known, the time interval between transmitted laser pulse and return is used to calculate target range and subsequently to determine the three-dimensional (3-D) location of the laser point measurement. Each of these geometric components has the potential to introduce error into the final elevation measurement. Further, complex topography and multiple reflections may combine with system geometry to induce additional measurement errors.
An emitted laser pulse diverges as it travels away from the source. Sensor beam divergence will typically vary from 0.2 to 0.3 mrad. A lidar beam diverges at a known rate. For example, a target spot radius of 0.2–0.3 m is the result at a range of 1000 m (Reference BaltsaviasBaltsavias, 1999a). Until recently, a typical lidar sensor provided only discrete returns, which represented peaks in the backscattered lidar illumination and its associated magnitude (intensity). Newer lidar systems have the ability to record discrete samples of the entire back-scattered illumination at significantly faster sampling rates (1–2 GHz). The resulting dataset is a discretely sampled waveform of the backscattered illumination (e.g. full-waveform lidar), which allows more accurate retrievals of surface elevations in a multiple reflection environment (Reference Mallet and BretarMallet and Bretar, 2009). In addition to multi-range measurements, physical properties of the targets included in the reflected illumination may be revealed by analysis of the shape of the sampled backscatter sequence (Fig. 3), and radiometric calibration allows for calibrated reflectance or amplitude retrieval (Reference WagnerWagner, 2010). The waveform or multiple returns allow for mapping of vegetation height, structure and/or understory in addition to enhancing the ability to map sub-canopy terrain.
The recent uptake in ground-based tripod lidar systems (terrestrial laser scanner, TLS) has significantly changed the acquisition of high-precision terrain data across scientific communities, including snow science (e.g. Reference ProkopProkop, 2008; Reference Prokop, Schirmer, Rub, Lehning and StockerProkop and others, 2008; Reference SchaffhausserSchaffhauser and others, 2008; Reference Mott, Schirmer and LehningMott and others, 2011; Reference Egli, Jonas, Grünewald, Schirmer and BurlandoEgli and others, 2012). Although these scanners are still expensive (>$150 000), a single scanner can provide a pool of researchers access to data year-round. TLS systems are less expensive to operate, and provide higher-accuracy (5 mm) and repeatable data with significantly fewer sources of error than airborne systems. Although a precision GPS is still required to achieve accurate point-cloud georegistration simultaneously or during post-processing, the use of an IMU and tracking motion associated with aircraft trajectories is eliminated, avoiding the largest source of measurement error in airborne lidar surveys. An inclusive scanner is tightly integrated with its scanning head mechanism and range encoder. TLS surveys do have their limitations. Scanners are line-of-sight and often range-limited (400–1000 m), which requires multiple scan location occupations to create complete coverage of the area of interest. For large survey areas, TLS surveys are typically less productive than airborne systems.
Photon-counting, or Geiger mode, lidar is seeing increasing use and development (Reference DegnanDegnan, 2002; Reference Harding, Shan and TothHarding, 2009). These systems require much lower return energy power and therefore can employ a low-power micropulse laser, with its attendant long laser life and inherent eye safety. A photon-counting system measures the time-of-flight of every detected photon while assuming that each photon originated from the lidar system. Much noise is recorded due to solar-generated photons, but the photon density histogram reveals a return energy profile similar to a waveform from a conventional pulsed system, even at solar noon through moderate cloud cover (Fig. 3; Reference Harding, Shan and TothHarding, 2009).
2.2. Lidar system parameters
2.2.1. Common airborne and ground-based system parameters
The spatial resolution of the point data is often quantified as the average point density per square meter or by an average ground-point spacing, though occasionally it is represented as the smallest elevation contour interval that can be mapped from the data. Factors influencing the spatial point density at ground level from airborne systems include the scan pattern, scan rate, swath width, pulse-repetition frequency (PRF), aircraft speed and target range or aircraft height (Reference BaltsaviasBaltsavias, 1999a; Reference RenslowRenslow, 2012; Fig. 2).
2.2.2. Airborne system parameters
PRF is the primary determinant of point spacing and on older lidar systems can be a limiting factor. Newer commercial systems can achieve pulse rates >400 kHz, allowing for very dense laser shot patterns (Reigl, 2011). The scan rate is the angular velocity of the oscillating mirror that directs the outgoing laser pulse (commonly 0–100 Hz) and combines with the PRF to determine across-track point spacing (Reference Wehr and LohrWehr and Lohr, 1999).
Several scan patterns are currently in use for airborne sensors, and more complex scan patterns are being developed for acquisition of specific point densities or spatial coverages. The most common are parallel or Z-shaped bidirectional scans. Palmer elliptical scans are also used. They provide fore and aft pointing angles in addition to the across-track directions achieved by the bidirectional scan patterns and therefore theoretically provide more opportunities for forest canopy penetration. Scan angle is the primary control on forest canopy penetration, with smaller scan angles (determined from zenith) resulting in less oblique shot angles and a higher probability of canopy penetration and ground detection (Reference Lovell, Jupp, Newnham, Coops and CulvenorLovell and others, 2005).
Along-track point spacing for airborne systems is controlled by the aircraft ground speed and the scan rate, with the maximum spacing occurring coincident with the edge of the field of view (Reference BaltsaviasBaltsavias, 1999a). Adjacent swaths are overlapped to provide additional point density along the swath margins (Reference Wehr and LohrWehr and Lohr, 1999), with swath overlaps of 20–30% recommended (e.g. Reference LatypovLatypov, 2002). In practice, ground speed and scan period are often constrained so that the along-track spacing is consistent with the across-track spacing.
The width of the scanned swath is of great importance for mission planning purposes. Wider swaths allow greater areal coverage with fewer flight strips and therefore can significantly decrease the data collection cost. Swath width depends primarily on the scan angle of the scanner system and the aircraft altitude. Increasing swath width via aircraft height comes at the expense of laser point spacing and/or range accuracy. At larger scan angles, the probability of forest canopy penetration decreases, with likely reductions in achievable ground-point spacing. Additionally, pointing accuracy decreases at the edge of the scan, with greater accuracy reductions at wide scan angles and high scan frequencies. This effect is magnified at higher flight altitudes, and scan angles are commonly decreased for high-altitude surveys.
2.2.3. Ground-based system parameters
TLS systems share many system parameters with airborne systems: pulse repetition frequency, scan rate, beam divergence and ground-point spacing are all the same or analogous. Ground-based scanners typically use a rotating or oscillating mirror to scan across either the vertical or horizontal scan angle (), analogous to the cross-track scan in airborne systems (Riegl: http://www.riegl.com/; Optech: http://optech.ca; Leica Geosystems, 2012: http://www. leica-geosystems.us; Fig. 2). TLS systems also typically rotate about the axis opposite the mirror axis with a rotation speed determined by the step angle increment and the time to complete a vertical scan line. Additionally, an integrated camera is common (either internal or a coupled external single-lens reflex), providing the capability to assign true color to the resulting point cloud and aid in point classification and scan registration (e.g. Reference Ullrich, Schwarz, Kager, Gruen and KahmenUllrich and others, 2003).
2.2.4. Post-processing
The collected data, after georegistration, are represented as (x,y,z) points, where z is a vector or matrix representing a single return, multiple returns or a full waveform, with corresponding attributes (e.g. return intensity or transmit pulse information). Sophisticated systems will geolocate the entire waveform or all returns, with the resulting z being a spatial vector (e.g. Reference Blair, Rabina and HoftonBlair and others, 1999). The raw data are commonly classified or filtered in order to group all points belonging to a surface of interest, i.e. top-of-vegetation (vegetation) or bare-earth (ground/snow surface) returns.
Most classification algorithms are proprietary to commercial software packages that are specific to lidar systems or laser-mapping contractors, which creates a potential source of error that is unavailable to the data user. However, in recent years there has been a surge in open-source software libraries tailored to lidar processing, analysis and format translation (e.g. LibLAS, PDAL and LAStools). Most data classification is done via automated algorithms that are monitored manually, often in combination with co-collected camera imagery (Reference Wehr and LohrWehr and Lohr, 1999).
Once point processing is complete, file formats are typically exported to a well-known format for larger distribution or further analysis or processing in other remote-sensing software packages. Point formats can typically be broken into two categories, vendor-specific and specifications and standards, though ASCII text files of (x,y,z) or (x,y,z,intensity) are commonly used for smaller datasets. Vendor formats are typically formats developed by laser manufacturers for the processing and delivery of data derived from a specific sensor (e.g. Riegl, Leica or Optech). All major lidar manufacturers produce their own proprietary format for use within their own processing software environment. For larger distribution within the geospatial and scientific communities there are a series of more generalized formats available for point data distribution. The LASer (LAS) specification released by the American Society for Photogrammetry and Remote Sensing (ASPRS) is the predominant method and de facto standard for distributing point data (http://www.asprs.org/a/society/committees/standards/LAS_1_4_r12.pdf). LAS is a binary file format designed for transmitting point-cloud data that is supported by all contemporary lidar processing software tools. The current version of the LAS specification is LAS 1.4 and is managed by the LAS working group under the oversight of the lidar division of ASPRS (Reference GrahamGraham, 2012). It should be noted that LAS is a file specification and does not adhere to the typical rigorous requirements of other well-known standards in the geospatial community, such as reference implementations, released software libraries and open format discussions. An alternate file format recently available is E57, which is a standard designed as a more comprehensive format for 3-D laser-imaging systems and was released by ASTM (http://www.astm.org/commit.E57.htm). E57 is a relatively new format and does not have as large an uptake as LAS. However, E57 was developed under a consortium of vendors and commercial entities and has a more adaptable position to be a comprehensive format in the lidar market.
The companion format to LAS is LAZ (LASzip), a lossless compression format based on the LAS specification. The LAZ format utilizes the open-source LASzip software library to reduce file sizes to ∼7–20% of the original LAS file size. LAZ is not an ASPRS-sanctioned format but rather a product developed independently based on the ASPRS LAS specification (Reference IsenburgIsenburg, 2013).
After classification of the point data, the snow-off and snow-on datasets must be co-registered to ensure proper alignment, commonly using building roofs or terrestrial features that are free of snow in both datasets (e.g. Reference Besl and McKayBesl and McKay, 1992; Reference Fowler and KadatskiyFowler and Kadatskiy, 2011). After co-registration the snow-free ground elevations can be subtracted from the snow-covered elevations to derive snow depth.
2.3. Snow-depth calculation
Calculation of snow depth involves subtraction of snow-free from snow-covered surface datasets, analogous to the techniques used for land surface change detection (Fig. 4). Specialized processing software packages can directly subtract two point-cloud datasets (e.g. Cloud Compare, http://www.danielgm.net/cc/), producing a difference cloud, which can be analyzed point-wise or gridded.
For users without specialized software and/or where gridded snow depth is the desired product, the most efficient method of subtracting the two datasets is to interpolate both snow-free and snow-covered point clouds to common grids and subtract the grid values. An alternate method is to convert the snow-free elevations to a grid dataset, extract the grid values below each snow surface elevation point measurement and perform a point-wise subtraction (Reference Deems, Fassnacht and ElderDeems and others, 2006). This method produces a point dataset that maintains the measurement locations and distribution of point spacings inherent in the snow-on point cloud. The creation of the snow-free grid dataset, or digital elevation model (DEM), involves interpolation and thus potentially induces some error, dependent on the terrain surface complexity and the point spacing. The high spatial resolution of the lidar point data commonly allows use of a simple interpolation scheme, such as inverse-distance weighting, with minimal introduced error. However, steep terrain or areas with substantial ground cover or forest understory may require more sophisticated interpolation schemes to minimize interpolation errors, which can be decimeter-scale in sloped or shrub-covered areas (Reference Bater and CoopsBater and Coops, 2009). Terrain with significant vertical displacements, such as cliff bands, presents challenges to DEM generation by simple interpolation; however, GIS techniques, such as barrier or break-line delineation, can be used as needed (e.g. Reference BrieseBriese, 2004).
For analyses involving grid generation, the grid element size should be of similar magnitude to the average point spacing of the filtered elevation point dataset in order to minimize scaling concerns and smoothing. In general, a small number of nearest neighbors should be used to interpolate each grid value in order to minimize smoothing errors. The highest degree of smoothing will occur in areas of lower point density, such as where heavy forest cover exists or potentially in areas where terrain dictates a higher flight elevation or when turbulence of the aircraft creates voids in the planned sampling pattern (Reference RenslowRenslow, 2012).
2.4. Integration with other sensors
High-resolution digital photography is often acquired concurrently with lidar data in both ground-based and airborne applications, and its use during the classification process can greatly improve data quality and reduce ambiguity of return classification in areas where automated filters have difficulty. If the relative geometry of the camera and lidar systems is known, orthorectified photographs can be used to provide color-mapping information for visualization of the lidar point cloud (Fig. 1). The snow-covered imagery will also elucidate snowdrift and scour features, enabling a qualitative assessment of the final snow-depth map.
In areas or seasons with partial snow coverage, snow-covered area products, such as fractional snow-cover maps from the NASA Moderate Resolution Imaging Spectroradiometer (MODIS) or the NPOESS Preparatory Project (NPP) Visible/Infrared Imaging Radiometer Suite (VIIRS) (e.g. Reference Painter, Rittger, McKenzie, Slaughter, Davis and DozierPainter and others, 2009), could be used for mission planning, allowing mapping of snow-free area to be minimized and thereby reducing flight time.
Reference Fowler and KadatskiyFowler and Kadatskiy (2011) have pioneered the integration of TLS survey data to rectify Airborne Laser Scanner (ALS) survey data. A common error assessment technique is to adjust the millions of points in an ALS survey to match a small number of surveyed ground-control points, which limits the potential final survey accuracy. By using a terrain model derived from a coincident TLS survey to perform a least-squares fit, the accuracy of the ALS-derived terrain model can be greatly improved (Reference Fowler and KadatskiyFowler and Kadatskiy, 2011). Recently, national programs and specifications have sought to address overall accuracy in airborne lidar collects, and a trend is appearing that requires more rigorous or enhanced validation methodologies for lidar collections (Reference HeidemannHeidemann, 2012).
3. Error Sources in Lidar Mapping
Here we briefly review error sources common to most laser-mapping applications. For more detail, the reader is referred to Reference BaltsaviasBaltsavias (1999a), Reference Wehr and LohrWehr and Lohr (1999), Reference Hodgson and BresnahanHodgson and Bresnahan (2004) and Reference GlennieGlennie (2007). Errors in snow-depth calculations from volumetric scattering of laser light by near-surface snow layers are discussed in Section 4.
Accuracy and resolution in lidar mapping are often viewed differently by data users and providers. Airborne lidar system accuracies are specified separately for the horizontal (x,y) and vertical (z) as the consequence of the above-ground view geometry is that vertical error is governed largely by errors in the range measurement, while horizontal error is primarily determined by attitude measurement errors (Reference GlennieGlennie, 2007). Vertical accuracy (often referred to as vertical tolerance by lidar contractors) refers to the degree to which the measurement reports the true elevation of the laser target. Horizontal accuracy describes the planimetric locational accuracy of each ground spot, typically reported to be on the order of 1/1000 of flight height (Reference Hodgson and BresnahanHodgson and Bresnahan, 2004). Horizontal accuracy varies across-track with pointing angle – reported values are commonly the average of maximum and minimum pointing angles (edge of swath and nadir; http://www.airborne1.com/ LiDARAccuracy.pdf). Additionally, accuracies are usually specified for flat open terrain of 10–20% reflectance and typically are 1σ specifications based on statistical sampling of measurements, meaning that 68% of measurements will fall within the specified error bounds.
Horizontal and vertical errors in TLS systems are strongly coupled, unlike ALS systems, obviating the need for separate error descriptions. Geolocation errors are the primary determinant of 3-D position error in TLS surveys.
Horizontal resolution refers to the nominal ground-point spacing between laser spot centers. Actual ground-point spacing in ALS surveys will be subject to the horizontal accuracy influenced by terrain geometry, but resolution is commonly a primary planning target specified to achieve specific project goals and serves as a primary driver for choosing flight and system parameters for lidar data collection. Ground-point spacing in TLS surveys is determined by the pulse and angular scan rates and increases with distance from the scan location.
3.1. Boresight calibration
For airborne or mobile systems, the boresight offset, the angular offset between the laser scanner and the IMU, must be calibrated. Boresight errors, as with all angular errors, increase with range to target and can contribute on the order of 3–5 cm to the vertical error and 10–20 cm to the horizontal error at 1000 m above ground level (Reference GlennieGlennie, 2007). Bore-sight calibration is usually conducted using overlapping lidar swaths, which are often flown in different directions and at different altitudes over a known target, such as the airport runway. Reference GlennieGlennie (2007) demonstrates that a least-squares boresight adjustment is superior to manual adjustment.
3.2. Global navigation satellite system (GNSS) and inertial measurement unit
The GNSS/GPS and IMU systems provide positional and platform orientation data from which the location and attitude of the laser sensor can be derived and thus the locations of the ranged ground elevations can be determined precisely (Reference Wehr and LohrWehr and Lohr, 1999). The GPS and inertial navigation system (INS) must be time-rectified with the laser scanner so that all laser range measurements are tied to the appropriate positional data. GPS accuracy enhancements, such as DGPS, real-time kinematic (RTK) or Precise Point Positioning (PPP), are required to achieve sub-decimeter positional accuracy (Reference Kaplan and HegartyKaplan and Hegarty, 2006). DGPS and RTK require a nearby ground reference GPS station, the distance to which affects the accuracy of the position solution (Reference KingKing, 2009). The time series of GPS locations is then corrected (either in real time for RTK or at a later time for DGPS or PPP techniques) using the time-series offsets recorded by the ground station. Range error magnitudes due to GPS/IMU are typically on the order of 1–2 dm (Reference BaltsaviasBaltsavias, 1999a; Reference MaasMaas, 2002; Reference KingKing, 2009).
3.3. Terrain-induced errors
Positional error due to terrain slope occurs via two mechanisms. First, errors in the horizontal (x,y) directions will induce an apparent error due to the uncertainty of the planimetric position of the measured elevation (Fig. 5a). This effect of vertical error dependence on horizontal accuracy can cause the measured point to appear above or below the actual terrain surface due to errors perpendicular to the contour direction. Horizontal error for airborne data is commonly estimated to be on the order of 1/1000 of aircraft flight altitude (above ground level) (Reference Hodgson and BresnahanHodgson and Bresnahan, 2004).
The second slope-induced error is due to the spreading of the laser spot on the inclined surface (Fig. 5b). This effect will spread the time distribution of the returned pulse, increasing the ‘rise time’ for the return to reach the intensity threshold for return signal registration and thus increasing the recorded range distance. For a 45° slope with a flight height of 1000 m, this ‘time-walk’ effect can induce a vertical error close to 50 cm (Reference BaltsaviasBaltsavias, 1999a). Smaller beam divergence angles produce smaller ground footprints, minimizing this effect. Flight planning can also help minimize time-walk by accounting for areas of steep terrain by assessing the laser-scanning and terrain geometry and orienting flight-lines to minimize the number of oblique-incident laser shots on steep slopes (Fig. 6). Errors induced due to scanner/terrain geometry combinations are not commonly assessed or quantified and therefore are usually not included in error budgets despite the potentially substantial contribution to the overall accuracy of the dataset, particularly in rough terrain (Reference Schaer, Skaloud, Landtwing and LegatSchaer and others, 2007). Incidence angle error contributions can be substantial in TLS scans as well despite their much smaller laser footprint and should be considered when choosing scan locations (Reference Soudarissanane, Lindenbergh, Menenti and TeunissenSoudarissanane and others, 2009). Fine-scale surface roughness can impart local slope angle variations at the resolution of TLS scan patterns, which will appear as increased millimeter- to centimeter-scale noise on the broader terrain surface.
3.4. Vegetation-induced errors
The ability of lidar to map both forest canopy and ground or snow surface elevations in one survey is one of the more attractive features of the technology. Accurate sub-canopy mapping is contingent on a sufficient number of laser shots reaching the ground and returning to the sensor directly. The number of successful ground hits and therefore the final surveyed point density decreases inversely with canopy cover density. Reference Reutebuch, McGaughey, Andersen and CarsonReutebuch and others (2003) compared lidar and manually surveyed elevations in several areas that are open, forested or forested with understory. Their results showed a relatively minor sub-canopy decrease in ground-point density, with point spacings on the order of 1 point m−2 in old growth Douglas Fir forest. The accuracy of the measured point elevations was degraded by an order of 10 cm by the forest or understory due to reduction in the strength of the return signal.
The specific relationship between ground return point density and canopy density is determined by the forest cover type (and therefore canopy and understory structures), laser spot size, the laser pulse rate and the scan angle of the laser sensor. The effects of canopy screening are minimized by increased pulse rates and decreased scan angles. Higher pulse rates on the order of 50–500 kHz provide more laser shots per square meter and thus an increased probability of successful canopy penetration. Smaller scan angles increase the probability of canopy penetration by reducing the number of individual trees that a single laser shot must penetrate (Fig. 1). Smaller laser spot sizes due to smaller beam divergence or lower flight altitude also increase the likelihood of canopy penetration and ground detection by increasing return energy. However, sub-canopy elevation mapping is of comparable accuracy to surveys conducted in open areas (e.g. Reference Reutebuch, McGaughey, Andersen and CarsonReutebuch and others, 2003; Reference HodgsonHodgson and others, 2005). In snow-covered landscapes, the reduction of understory vegetation due to snow burial combined with the high reflectivity of the snow surface at lidar wavelengths of 1064 nm or shorter should allow for increased sub-canopy ground-point density and accuracy compared with snow-free collections. Ground-based systems operating at 1550 nm wavelengths will still benefit from reduced understory vegetation during the snow season, though the snow reflectance at that wavelength offers no advantage over bare ground.
3.5. The role of flight and survey planning in error control
Flight planning is critical to mission success, minimizing cost and data accuracy, especially in rough terrain (cf. Reference SaylamSaylam, 2009, for a brief overview of flight planning considerations). Terrain geometry, such as slope magnitude and aspect relative to the flight-line, interacts with the laser scan angle to affect the laser angle of incidence at a given ground surface location. In areas with very steep terrain (e.g. near-vertical to overhanging), flight tracks along or above the steep terrain features will result in terrain shadowing and gaps in laser returns, with important ramifications for any subsequent interpolation or gridding of the surface elevations (e.g. Reference Hopkinson and DemuthHopkinson and Demuth, 2006). Proper flight planning will minimize terrain shadows and the number of points collected with poor geometry.
Figure 6 illustrates the influence on laser shot geometry of flight-line orientation relative to terrain. A flight-line oriented along a ridgeline will produce a preponderance of laser shots with oblique incidence angles (flight-line 1), while a flight-line offset from the ridgeline will provide a much higher proportion of laser shots with acute incidence angles, minimizing errors due to beam spreading and forward scattering.
A similar terrain geometry issue pertains to TLS surveys. Scan locations that produce predominately downslope (high incidence angle) laser geometries will have greater errors than those that achieve incidence angles closer to slope normal (Reference Soudarissanane, Lindenbergh, Menenti and TeunissenSoudarissanane and others, 2009). In addition to incidence angle geometries, TLS survey plans should explicitly allow for one or more tie points (usually reflectors) and/or planar features to be common between adjacent scans. Tie points are essential in establishing the registration of multiple scan positions and the georegistration of the data (Reference JacobsJacobs, 2005). Fewer tie points are required if scan positions and reflectors are located with high-precision GPS; with lower-quality GPS position information, a greater number of common tie points will be needed to enable consistent georegistration and mutual adjustment of overlapping scans through techniques such as multi-station adjustment (e.g. Reference Ullrich, Schwarz, Kager, Gruen and KahmenUllrich and others, 2003; Riegl, 2012).
3.6. Post-processing errors
After data collection and georegistration, the ‘point cloud’ of raw lidar return points is classified as ‘terrain’ or ‘non-terrain’returns. The point-cloud classification is often highly automated. The classification can be accomplished using any number of (often proprietary) algorithms, but commonly involves segregating terrain points through an iterative process that evaluates the deviation of individual points from a surface generated from nearby points. Classification thresholds based on these deviations can be somewhat subjective, requiring a priori knowledge of the survey area and/or manual supervision using ancillary data such as digital orthophotographs or existing lower-resolution DEMs (Reference Zhang, Chen, Whitman, Shyu, Yan and ZhangZhang and others, 2003). Useful overviews of classification methods are provided by Reference LiuLiu (2008), Reference Evans, Hudak, Faux and SmithEvans and others (2009), Reference Meng, Currit and ZhaoMeng and others (2010) and Reference TinkhamTinkham and others (2011).
Misclassification of points can induce errors in the final elevation surface. Because the error magnitude depends on many factors (e.g. the type of filter used, the accuracy of the measured elevation, the elevation of the vegetation above the terrain surface and the terrain geometry), the contribution of classification errors to the overall surface accuracy will vary widely. It is clear, however, that successful application of classification algorithms is critical to the accuracy of the final elevation surface or snow-depth map (Reference Evans and HudakEvans and Hudak, 2007). Though most classification routines can be configured to run with minimal operator interaction, manual inspection of the geolocated point cloud and of subsequent filtered products in conjunction with collocated imagery can be important for qualitative validation of automated routines.
4. Snow Surface Interactions
Snow is composed of particulate ice, air, liquid water and impurities such as dust, soot and organics. The snow particle size is described in terms of optical grain radius or specific surface area (SSA; m2 kg−1); here we will use the optical grain radius. Snow optical grain radius usually ranges from ∼50 to >1000 μm (Reference Painter, Dozier, Roberts, Davis and GreenPainter and others, 2003). These particles create a volume-scattering medium with penetration of radiation into the snowpack, which depends primarily on wavelength and particle size.
Most airborne lidar acquisitions are performed with laser systems that operate at wavelength λ = 1064 or 532 nm, while many TLS systems operate at λ =1550 nm (Reference BaltsaviasBaltsavias, 1999b). In the near-infrared (NIR) wavelengths (1064 nm lidar), ice is moderately absorptive and therefore the snow reflectance is most sensitive to grain size, so penetration is less than a few centimeters (Reference WarrenWarren, 1982; Fig. 7). In the visible wavelengths (e.g. 532 nm lidar), absorption is much lower, resulting in higher snow reflectance but greater penetration of the incident laser light into the snowpack. In the shortwave infrared (e.g. 1550 nm lidar), absorption by ice is much stronger and the snow reflectance is <10%.
The optical properties of ice are described by the complex refractive index, m:
where n, the real part, describes refraction and k, the imaginary part, describes absorption (also termed the absorption coefficient). Table 1 shows values of k for ice at common lidar system wavelengths (from Reference WarrenWarren, 1984; Reference Warren, Brandt and GrenfellWarren and others, 2006). An understanding of the sensitivity of the light transmission and absorption at common lidar wavelengths to snow grain size and snow liquid water content is important for an assessment of the contribution of laser interaction with the snow surface to the uncertainties in snow-depth retrievals.
4.1. Light penetration of the snow surface
Measurements of the transmission of solar radiation through the near-surface layers of the snowpack, particularly at the wavelengths that interest us here, are extremely difficult if not impossible to make (see Reference Giddings and LaChapelleGiddings and LaChapelle, 1961). Therefore, the literature on these measurements and modeling of transmission of solar radiation through snow is sparse and our knowledge poorly constrained.
Reference Beaglehole, Ramanathan and RumbergBeaglehole and others (1998) and Reference PerovichPerovich (2007) measured the spectral transmission of solar radiation in snow in wavelengths from 350 to 900 nm in Antarctic and temperate snow covers. Both found reductions in transmission of >80% at all wavelengths by a depth of 2 cm in the snow columns, i.e. >80% of the signal in those wavelengths is reflected by the snow volume before reaching 2 cm depth. At 500 nm, wavelength transmission remained above 5% by 10 cm depth, whereas at 900 nm wavelength transmission of radiation at 4 cm depth dropped to ∼0.006. Given that k is ∼4.1 × 10−7 at this wavelength, whereas at 1064 nm it is order 10−6, we can consider that the same attenuation of radiation (scattered + absorbed) comes at a shallower depth (order 1 cm). We are unaware of any transmission measurements in snow at 1550 nm. Nevertheless, the strong absorption at 1550 nm will bound the transmission to less than that at 1064 nm. Therefore, the vast majority of reflected signal at 532 nm comes from the top 10 cm, at 1064 nm the majority comes from the top 1 cm and at 1550 nm likewise from the top 1 cm.
Light-absorbing impurities in snow generally consist of dust, carbonaceous particles and pollen, each of which has less forward single scattering than ice grains (Reference Wiscombe and WarrenWiscombe and Warren, 1980; Reference WarrenWarren, 1982). With lower forward scattering, snow with light-absorbing impurities should have less penetration into the snowpack than clean snow and therefore buffer the returned signal closer to the snow– atmosphere interface.
At 1064 nm, the light-absorbing impurities will have little effect on the spectral reflectance and intensity of returned signal given that the contrast in k for ice and impurities is relatively small and the proportion of ice is vastly greater than most impurity concentrations (Reference WarrenWarren, 1982). Only when impurity concentrations reach the higher levels of parts per thousand does the spectral reflectance begin to drop (Reference Singh, Kulkarni and ChaudharySingh and others, 2010; Reference PainterPainter, 2011). Light-absorbing impurities have their most powerful impact on reflectance and returned signal for 532 nm lidars because the contrast in k between ice and the impurities is generally greatest in the visible wavelengths.
Liquid water content should have a similar effect to a coarsening of grain size, whereby the optical depth of the snowpack decreases (increasing transmission) but the absorbing path length increases (decreasing transmission). Therefore, given the bounding measurements, we estimate that 97% of attenuation occurs in the top 1 cm of the snowpack. The presence of liquid water at and near the snow surface, though imparting minimal range error, will strongly decrease return energy for NIR wavelengths, potentially reducing point densities due to dropouts and limiting range for TLS systems.
In summary, the sensitivity of the error budget to scattering depth is less than the potential contributions of scattering from fine-scale surface roughness features, temporal location mismatches and horizontal geolocation errors.
4.2. Anisotropic reflectance
Snow is strongly forward-scattering and the proportion of forward scattering increases with wavelength (as does absorption) (Reference WarrenWarren, 1982; Reference Painter and DozierPainter and Dozier, 2004). Likewise, the anisotropy of the angular distribution of reflected radiation increases to the forward direction as snow grains increase in size (Fig. 8). These factors combine to reduce lidar return strength for laser shots oblique to terrain slope and for snowpacks of increasing grain size, thus increasing the rise time for return registration and increasing the likelihood of dropped/missed returns.
For example, for a downslope laser shot, a fine-grained winter snowpack at λ = 532 nm will have the strongest return signal, while a coarse-grained wet spring snowpack at λ = 1550 nm will have the weakest return signal. The decrease in return signal strength due to forward scattering induces the same time-walk effect as slope-induced beam spreading, slowing the rise time for the return energy to reach the detection threshold for discrete-return systems and contributing to the vertical error budget.
5. Recommendations for Lidar Snow Surveys in MOuntainous Areas
5.1. Airborne surveys
Much of the potential for error exists in the geometric relationships between the aircraft dynamics, the scanning system and the irregular ground surface. Careful flight planning is critical for minimizing these errors. Most lidar contractors use flight planning models to constrain many of the laser geometry parameters, but it is rare that terrain considerations other than flight altitude (above ground level) are considered – nominal ground-point spacing is usually the primary planning metric. We recommend that investigators seeking lidar acquisitions be proactive and involved in the planning process, especially if the target area contains steep terrain.
Of particular concern in high-relief areas is the constraint imposed on pulse repetition frequency (PRF) by aircraft altitude. If the time of flight of an individual laser pulse is longer than the interval between pulses (PRF), the pulse echo will arrive after the next pulse is transmitted, complicating allocation of the return pulse to a specific transmission. This multiple-time-around (MTA) or multiple pulse in air (MPiA) effect requires lower PRF at higher flight altitudes, unless the lidar system is specifically designed to resolve range ambiguities (Reference LemmensLemmens, 2011).
The maximum unambiguous measurement range is the time of flight allowed between laser pulses (Reference Rieger and UllrichRieger and Ullrich, 2011):
For example, a lidar system operated at 400 kHz PRF has a maximum unambiguous measurement range of 375 m, far less than the maximum measurement range and obviously a severe flight altitude restriction in mountain basins where relief can be much greater.
There are a number of methods for resolving range ambiguities with pulsed lidars. Using two scanners with different wavelengths or pointed at different halves of the ground swath, the PRF can be effectively doubled. As a post-processing technique, an existing DEM can be used to establish expected ranges, but this technique is vulnerable to errors in the existing terrain model and application is limited to areas with existing DEMs (e.g. Reference Roth and ThompsonRoth and Thompson, 2008). More sophisticated methods involve real-time or near-real-time statistical processing to establish correlations between transmitted and received pulses (e.g. Reference Hiskett, Parry, McCarthy and BullerHiskett and others, 2008) or to minimize noise within a single scan line induced by an imposed minor variation (or ‘jitter’) in PRF (Reference Rieger and UllrichRieger and Ullrich, 2011). Future lidar systems may employ pulse-compression or pulse-tagging techniques for unambiguous return assignment. In high-relief terrain, without using a system designed to resolve range ambiguities, reduced ground-point spacing (due to reduced PRF) must be acceptable or flight planning must incorporate multiple flight altitudes to constrain measurement ranges to less than the maximum unambiguous range.
For snow-depth mapping in mountainous forested terrain, high laser pulse rates (>50 kHz) and scan angles up to order 30° (15° each side) are desirable to minimize slope-induced errors and to maximize canopy penetration and thereby maximize ground/snow surface point densities. Flight-lines require careful planning so that laser shot/slope angle geometries remain favorable, minimizing oblique down-slope shots (Fig. 6). In very complex terrain, significant swath overlap may be required to ensure that the planned point spacing is achieved. Flight planning requirements for complex mountain terrain, especially the potential for complicated flight-lines, are likely to significantly increase flight time and data acquisition costs over survey criteria that are adequate for flat unvegetated sites.
5.2. Ground-based surveys
Typical TLS systems that are optimized for topographic scanning at moderate- (100–500 m) to long-range (500– 1500 m) scales operate at 1550 μm in the NIR and are eye-safe. When deployed to characterize snow and ice they are significantly range-limited (<150 m) due to the very low reflectance of snow and ice at 1550 μm, especially if the snow surface is wet or consists of large grains (cf. Reference ProkopProkop, 2008; Fig. 7). An alternative is to utilize a 1064 μm wavelength for increased range and reflectivity, although utilization of a 1064 μm sensor is limited to a select few commercially available sensors (e.g. Optech ILRIS-LR) and some custom sensors (e.g. Riegl VZ-6000). TLS systems with a 1064 μm wavelength are not eye-safe at close range using common pulse and scan rates.
Georegistration of on-snow scanning using TLS is fundamental to successful acquisition. Achieving a stable tripod set-up can be difficult on a snow surface, especially during the melt season or in areas with non-supportable snow. In a typical TLS operation, scans are registered to one another utilizing in-field reflectors and precise GPS coordinates. Establishing stable reflector control points can be difficult due to snow settlement or deformation by melt, creep, metamorphism or solar influence. An alternative approach is to utilize a single scan position with an RTK GPS receiver on a TLS system with internal compass abilities (e.g. 2012 Riegl VZ-series). Utilization of an internal compass and on-board RTK reduces the need for destruction of the snowpack for establishing control and allows for single scan position scanning and post-processing registration for precise scan alignment.
6. Conclusions
Lidar surveying is used increasingly as a data source for high-resolution terrain data. Repeat surveys including one snow-free collection and one or more snow-covered collections allow the calculation of snow depth over sizeable geographic areas at fine resolutions and with high vertical accuracy, giving numerous advantages over manual survey or larger-footprint sensors. Lidar can map snow depth, the primary source of spatial variability in SWE, at high spatial resolutions. It therefore has tremendous potential for snow hydrology and avalanche science applications and is sure to have a significant impact on snow science in the years to come.
Acknowledgements
Part of this work was supported by the National Snow and Ice Data Center, the national Oceanic and Atmospheric Administration (NOAA) Western Water Assessment and a Cooperative Institute for Research in Environmental Sciences (CIRES) Innovative Research Project grant. Part of this work was performed at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Part of this work was performed at the US Army Corps of Engineers, Cold Regions Research and Engineering Laboratory, Hanover, NH.