[go: up one dir, main page]
Next Article in Journal
Compressive Video Recovery Using Block Match Multi-Frame Motion Estimation Based on Single Pixel Cameras
Previous Article in Journal
SHM-Based Probabilistic Fatigue Life Prediction for Bridges Based on FE Model Updating
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Sensors
Volume 16
Issue 3
10.3390/s16030313
sensors-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Radiometric Calibration of a Dual-Wavelength, Full-Waveform Terrestrial Lidar

by
Zhan Li
1,2,*,
David L. B. Jupp
3,
Alan H. Strahler
1,*,
Crystal B. Schaaf
2,
Glenn Howe
4,
Kuravi Hewawasam
4,
Ewan S. Douglas
5,
Supriya Chakrabarti
4,
Timothy A. Cook
4,
Ian Paynter
2,
Edward J. Saenz
2 and
Michael Schaefer
3,6
1
Department of Earth and Environment, Boston University, 675 Commonwealth Avenue, Boston, MA 02215, USA
2
School for the Environment, University of Massachusetts Boston, 100 Morrissey Blvd., Boston, MA 02125, USA
3
CSIRO Land & Water, GPO Box 1666, Canberra, ACT 2601, Australia
4
Department of Physics and Applied Physics, University of Massachusetts Lowell, 600 Suffolk Street, Lowell, MA 01854, USA
5
Department of Astronomy, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA
6
Precision Agriculture Research Group, School of Science and Technology, University of New England, Armidale, NSW 2351, Australia
*
Authors to whom correspondence should be addressed.
Sensors 2016, 16(3), 313; https://doi.org/10.3390/s16030313
Submission received: 1 December 2015 / Revised: 23 February 2016 / Accepted: 23 February 2016 / Published: 2 March 2016
(This article belongs to the Section Remote Sensors)
Browse Figures
Figure 1
<p>DWEL system response pulse. The two pulse peaks at the two wavelengths are aligned.</p> "> Figure 2
<p>Example of saturated pulse and saturation correction.</p> "> Figure 3
<p>Calibration data collection. (<b>a</b>) Three panels for calibration, from top to bottom: painted light gray panel, white Spectralon panel and painted dark gray panel; (<b>b</b>) DWEL was set up in stationary mode with the laser pointing along the measuring tape laid out on the floor; panels were placed at 33 range locations along the tape; (<b>c</b>) the green laser was used to point the infrared lasers to the center of the panels.</p> "> Figure 4
<p>Estimation and validation of calibration of 1064 nm data. In both rows, the left column (<b>a,c</b>) shows the calibration function as fitted to training data, and the right column (<b>b,d</b>) shows the fit to the validation data. First row (a,b): measured and modeled intensity normalized by reflectance. Second row (c,d): scatter plots of measured against modeled intensity. The vertical error bars in (a) and horizontal error bars in (c) are one standard deviation of measured intensities normalized by reflectance.</p> "> Figure 5
<p>Estimation and validation of calibration of 1548 nm data. In both rows, the left column (<b>a,c</b>) shows the calibration function as fitted to training data, and the right column (<b>b,d</b>) shows the fit to the validation data. First row (a,b): measured and modeled intensity normalized by reflectance. Second row (c,d): scatter plots of measured against modeled intensity. The vertical error bars in (a) and horizontal error bars in (c) are one standard deviation of measured intensities normalized by reflectance.</p> "> Figure 6
<p>Telescope efficiency <math display="inline"> <semantics> <mrow> <mi>K</mi> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> </mrow> </semantics> </math> of the two wavelengths.</p> "> Figure 7
<p>Errors in apparent reflectance. The first row (<b>a</b>–<b>c</b>) shows 1064 nm, and the second row (<b>d</b>–<b>f</b>) shows 1548 nm. The left column (a,d) is the deviation from calibration fitting with range. The middle (b,e) column is the deviation of validation points with range. The right column (c,f) is the histogram of deviations.</p> "> Figure 8
<p>Sensitivity of the <math display="inline"> <semantics> <mrow> <msub> <mi>ρ</mi> <mrow> <mi>a</mi> <mi>p</mi> <mi>p</mi> </mrow> </msub> </mrow> </semantics> </math> estimate on errors in return intensity and range. The image color shows relative error in the <math display="inline"> <semantics> <mrow> <msub> <mi>ρ</mi> <mrow> <mi>a</mi> <mi>p</mi> <mi>p</mi> </mrow> </msub> </mrow> </semantics> </math> estimate (estimate − measurement). The color map scale is unified for all images for comparison purposes, but the actual error ranges of the four images are different and given here: (<b>a</b>) <math display="inline"> <semantics> <mrow> <mi>δ</mi> <msub> <mi>ρ</mi> <mi>α</mi> </msub> </mrow> </semantics> </math> at 1064 nm [−0.928, 0.928]; (<b>b</b>) <math display="inline"> <semantics> <mrow> <mi>δ</mi> <msub> <mi>ρ</mi> <mi>r</mi> </msub> </mrow> </semantics> </math> at 1064 nm, [−0.226, 0.290]; (<b>c</b>) <math display="inline"> <semantics> <mrow> <mi>δ</mi> <msub> <mi>ρ</mi> <mi>α</mi> </msub> </mrow> </semantics> </math> at 1548 nm, [−0.574, 0.574]; (<b>d</b>) <math display="inline"> <semantics> <mrow> <mi>δ</mi> <msub> <mi>ρ</mi> <mi>r</mi> </msub> </mrow> </semantics> </math> at 1548, [−0.133, 0.154].</p> ">
Versions Notes

Abstract

:
Radiometric calibration of the Dual-Wavelength Echidna® Lidar (DWEL), a full-waveform terrestrial laser scanner with two simultaneously-pulsing infrared lasers at 1064 nm and 1548 nm, provides accurate dual-wavelength apparent reflectance (ρapp), a physically-defined value that is related to the radiative and structural characteristics of scanned targets and independent of range and instrument optics and electronics. The errors of ρapp are 8.1% for 1064 nm and 6.4% for 1548 nm. A sensitivity analysis shows that ρapp error is dominated by range errors at near ranges, but by lidar intensity errors at far ranges. Our semi-empirical model for radiometric calibration combines a generalized logistic function to explicitly model telescopic effects due to defocusing of return signals at near range with a negative exponential function to model the fall-off of return intensity with range. Accurate values of ρapp from the radiometric calibration improve the quantification of vegetation structure, facilitate the comparison and coupling of lidar datasets from different instruments, campaigns or wavelengths and advance the utilization of bi- and multi-spectral information added to 3D scans by novel spectral lidars.

1. Introduction

Light detection and ranging (lidar) is an active remote sensing technique using an instrument that emits coherent laser light. Targets along the laser transmission path scatter the light, and the lidar instrument records the travel time and intensity of the scattered light received by its detector. A new and important application of lidar is the quantification of forest structure, principally measures of the physical dimensions of trees, the amount and location of leaves and gaps between and within tree canopies, through the 3D information acquired by lidar instruments on different remote sensing platforms (terrestrial, airborne and spaceborne). Studies have shown the capability of lidars to facilitate our understanding and management of forest ecosystems [1,2,3,4] by describing the heterogeneity of forest structure in 3D space and its relation to forest function [5,6].
A lidar instrument measures the range to a target through the product of the speed of light and the one-way travel time of light between the instrument and the target. The travel time can be measured using either pulse ranging or continuous wave ranging techniques [1]. For lidar remote sensing of vegetation, which is the application of concern in this paper, pulse ranging lidar is much more commonly used [1], and thus, we will confine our discussion here to pulsing lidar, unless otherwise noted.
Of the spatial location and intensity (while the term “intensity,” as defined in optical physics, refers to the energy flow rate in W·sr−1 from a point source of emission, we will use “intensity” here to refer to a measure, normally in digital counts, of the response of the detector-amplifier-digitizer system to the return of energy from a laser pulse or power of a continuous wave scattered by a target into the aperture of the telescope of the lidar instrument and reaching the detector system), the two primary attributes of scattering events recorded by lidar, location by range and zenith and azimuth angles (or as resolved into Cartesian coordinates) has found wide use in the retrieval of forest structural parameters [7], such as diameter at breast height (DBH) [8], tree and canopy height [9,10,11], timber volume [12,13], Leaf Area Index (LAI) [14,15,16,17,18] and others [19,20].
More complete inference of vegetation structure requires using intensity, the other attribute of the return signal. However, intensity information does not provide straightforward interpretation and has been underutilized. Intensities in digital counts output by lidar instruments neither give actual backscattered energy from targets nor relate directly to target physical properties. Accordingly, they are usually processed to remove electronic effects and normalize the decrease of observed intensity with range. The processed intensities thus provide the relative distribution of target return energy from which forest structural parameters can be inferred directly or through empirical regressions, such as gap fraction/probability [21], canopy height profile [22], basal area and aboveground biomass [23] and others [20]. Normalized intensity has also been used in target classification/recognition [14,24,25,26] and estimation of biochemical properties of vegetation [27,28].
Although these studies have documented the usefulness of lidar intensity values, their simple and empirical normalization is primarily arbitrary and only provides “relative reflectance” or “pseudo-reflectance” values through scaling of lidar intensities to adjust the contrast and overall “brightness” of lidar scans [29]. Empirical normalization limits the intercomparison of instruments, makes merging data from two or more instruments or scanning campaigns difficult and causes trouble in the interpretation of normalized intensity values. We need rigorous radiometric calibration to provide a value that is physically interpretable with regard to the properties of targets, such as surface reflectance.
The need for a consistent definition of calibrated lidar intensity is also driven by the recent design or fabrication of bi- or multi-spectral lidar instruments using lasers at different wavelengths or white lasers to exploit the spectral signatures of targets [30,31,32,33,34,35,36,37]. For example, the Dual-Wavelength Echidna Lidar (DWEL) instrument, which is the focus of this paper, uses two coaxial lasers at 1064 nm and 1548 nm wavelengths to differentiate leaves from branches, trunks and ground by taking advantage of the distinctive spectral response of leaves at the two wavelengths [30,38].
For airborne lidar scanning (ALS), recent studies have reviewed the physical concepts of return intensity and radiometric calibration [39,40,41]. Various calibration targets and procedures for ALS data have been proposed and evaluated for several calibration scenarios, including different scanning campaigns with the same instrument [42,43,44], different instruments at the same wavelength [45] or different wavelengths [46].
In terrestrial lidar scanning (TLS), calibrated intensities have been used to measure canopy structure. Examples include retrieval of the multiangle gap probability and then LAI [17], the clumping index estimate [18] and target classification with calibrated intensity alone [47] or along with pulse width from full-waveform data [48]. However, with the exception of a few recent studies on both pulse-ranging TLS [49,50,51] and continuous-wave ranging TLS [52], radiometric calibrations of TLS data are currently scattered among various application studies and are poorly documented or rely on undocumented proprietary calibration algorithms from instrument manufacturers. The radiometric calibration of TLS data faces unique challenges, including: (1) a very large variation in intensity with range that can induce saturation of the detector system by bright targets in the near field and reduced intensities that merge with the noise in the far range; and (2) strong telescopic effects, with defocusing that produces weak signals at near range.
This paper addresses these and other challenges for a dual-wavelength, full-waveform terrestrial laser scanner, the DWEL. No evaluation of the radiometric calibration of TLS data similar to the DWEL has yet been documented for different wavelengths. We present a simultaneous calibration of returns from DWEL’s two lasers, which demonstrates how calibration can ensure both radiometric and spectral fidelity in a unified process, thus providing a pathway for the calibration of other dual- and multi-wavelength terrestrial lidars that may now be in various stages of development and application.
In the paper, we begin with the theoretical lidar equation for canopy structure study using TLS data and then describe the processing of the DWEL waveform data. Our calibration model, based on a generalized logistic function for telescope efficiency and an inverse power fall-off with range, is fitted to stationary scans of panels with different reflectance values at different ranges. We conclude by evaluating the calibration accuracy of dual-wavelength point clouds from DWEL, as well as the sensitivity of the calibration accuracy to errors in both range and intensity measurements.

2. Physical Background

2.1. Basic Lidar Equation for Forest Canopies

The goal of our radiometric calibration is to obtain range-independent, instrument-independent and physically well-defined measurements for canopy structure modeling and estimation from returned power as detected and recorded by the lidar instrument’s optical and electronic systems. Previous studies [53,54] formulated lidar equations as a function of canopy structure parameters to model large-footprint airborne lidar waveforms, but did not identify a realizable quantity for lidar radiometric calibration. To establish the basis of the lidar calibration for the canopy structure study, we formulate the lidar equation from the basic scattering lidar equation [55] using Ross’s framework of radiation regime modeling of the vegetation canopy [56]. This lidar equation describes the interaction of the laser beam with vegetative elements and identifies the objective variable of the calibration (“apparent reflectance”).
Consider an angular voxel, an elemental volume enclosed by a laser beam between range r and r + Δ r from the lidar instrument. Vegetative elements inside one such angular voxel are modeled as a turbid medium composed of tiny thin facets of different orientations in space. We shall not specify the size and thickness of these facets nor their location inside the angular voxel [56]. Laser radiation incident to this angular voxel can be absorbed, reflected back toward the lidar instrument or transmitted through it without interaction with vegetative facets. Simulations of lidar waveforms with Monte Carlo ray tracing have shown that multiple scattering by vegetative elements in the canopy largely has no effect on return waveform shapes and contributes little to return energy, especially for the small laser beam divergence that we use here [57]. Thus, it is reasonable to assume only single scattering in the interaction between laser beams and vegetative facets.
The probability that a laser beam in a given direction reaches an angular voxel at range r without interaction with vegetative facets is given by the gap probability P g a p ( r ) :
P g a p ( r ) = e − G ∫ 0 r u L ( r ′ ) d r ′
where G is the Ross G-function, which describes the projection of a unit vegetative area in a given direction; u L ( r ) ( m 2 ⋅ m − 3 ) is the total upper side surface area of all tiny facets (no mutual-shading between facets, i.e., no clumping is considered) within a unit volume at range r along the laser beam [56]. Note that in the expressions below, we will consider only a single laser beam and omit the laser beam direction in P g a p ( r ) and G .
Let J 0 be the total outgoing laser radiation energy (units: J ) within an infinitesimal time, i.e., an impulse laser energy. The return energy Δ J (units: J ) from the angular voxel between range r to r + Δ r received by the telescope of the lidar instrument is:
Δ J = J 0 P g a p ( r ) ⋅ Δ β ⋅ K ( r ) ⋅ η s y s η a t m Δ β = ∫ Ω T = A T / r 2 1 π Γ ( r , Ω i → Ω v ) d Ω v ⋅ ∫ r r + Δ r u L ( r ′ ) d r ′
In these equations, J 0 P g a p ( r ) ( J ) is the laser radiant energy that reaches the angular voxel at range r ; Δ β (dimensionless) is the effective backscatter ratio of the angular voxel, i.e., the proportion of the incident laser radiation energy that is scattered back from the angular voxel into the solid angle subtended by the telescope receiving area; the expression   ( 1 / π ) Γ ( r , Ω i → Ω v ) ( s r − 1 ) is the area scattering phase function at range r , i.e., the portion of the radiant energy onto a unit area of vegetative facets in direction Ω i that is scattered in the direction Ω v within a unit solid angle, a term that contains both the radiative and structural characteristics of the vegetative facets [56]; Ω T ( s r ) is the solid angle subtended by the area of the telescope aperture ( A T ) of the lidar instrument from range r ; K ( r ) (dimensionless) is telescope efficiency at range r (see Section 2.1.3); and η s y s and η a t m (dimensionless) are transmission factors that account for energy loss due to the sensor system and atmosphere, respectively.

2.1.1. Apparent Reflectance

In the simplest case: (1) only one angular voxel along a laser beam is filled with vegetative elements at range r , that is P g a p ( r ) = 1 , and all laser radiation energy of this beam falls onto this voxel; (2) all of the vegetative facets in the angular voxel are Lambertian with the same diffuse reflectance ρ d (dimensionless), and all faces are orthogonal to the laser beam, that is ( 1 / π ) Γ ( Ω i → Ω v ) = ρ d / π ; and (3) all of the vegetative facets together fill the whole laser beam, that is ∫ r r + Δ r u L ( r ′ ) d r ′ = 1. The return energy from such a voxel is then:
Δ J = J 0 ⋅ ρ d ⋅ A T π r 2 ⋅ K ( r ) ⋅ η s y s η a t m
This simplest case is equivalent to an extended Lambertian panel (by “extended panel,” we describe a panel that is large enough to intercept the whole cross-section of a laser beam at a given range) with reflectance ρ d that fills the laser beam orthogonally at range r . The expression A T / ( π r 2 ) in Equation (3) is the proportion of the total hemispherical backscattering that is intercepted by the telescope aperture given this simplest case. Assuming η s y s and η a t m are constant with the range of our instrument operation and the laser wavelength in consideration, we may simplify Equation (3) by combining all constants into Φ 0 ( J ⋅ m 2 ), i.e.,
Δ J = Φ 0 ⋅ ρ d ⋅ K ( r ) r 2 Φ 0 = J 0 A T π η s y s η a t m
The reflectance ρ d , describing the physical properties of targets in the simplest case, is computable from the received energy and also is range independent and instrument independent. For the general case observed in reality, we cannot retrieve the anisotropic reflectance value of each vegetation facet, but may identify a variable called apparent reflectance, ρ a p p (dimensionless), to represent the overall structural and radiative characteristics of all of the vegetative facets as a whole illuminated by a laser beam. The ρ a p p for the general case is calculated in the same way as ρ d for the simplest case, i.e.,
Δ J = Φ 0 ⋅ ρ a p p   ⋅ K ( r ) r 2
Comparing Equation (5) to Equation (2), we also have,
ρ a p p = P g a p ( r ) Δ β A T / ( π r 2 )
Apparent reflectance, our range-independent and instrument-independent quantity, provides the objective variable for our lidar calibration for the canopy structure study for two reasons. First, ρ a p p can be modeled to derive canopy radiative and structural information. According to Equation (6), ρ a p p includes P g a p ( r ) , as determined by the structural characteristics of vegetative facets in the canopy, and Δ β , as determined by both the radiative and structural characteristics of the vegetative facets of the canopy. Second, ρ a p p can be interpreted as the reflectance value of a diffusely-reflecting, partly-absorbing panel filling the laser beam orthogonally that would return the same laser energy received by the lidar as the actual target at the same range. It can be realized as the ratio of (dark current corrected) Δ J , the lidar intensity from a target to Δ J w , the intensity from a white Lambertian panel (orthogonal to the laser beam, reflectance ρ w = 1 ) at the same range, as shown in Equation (7):
ρ a p p   = Δ J Δ J w ⋅ ρ w = Δ J Δ J w
Thus, ρ a p p is theoretically useful for the canopy structure study, as well as practically realizable for lidar calibration. It was introduced in Parkin et al. [58] and used in data interpretation by Jupp et al. [17].

2.1.2. Physical Interpretation of Apparent Reflectance

When the thickness of the angular voxel Δ r → 0 , we have the differential form of Equation (2) as:
J ( r ) = ∂ J ∂ r = J 0 P g a p ( r ) ⋅ β ( r ) ⋅ K ( r ) ⋅ η s y s η a t m β ( r ) = ∂ β ∂ r = ∫ Ω T 1 π Γ ( r , Ω i → Ω v ) d Ω v ⋅ u L ( r )
where J ( r ) ( J ⋅ m − 1 ) is the received laser energy from range r per unit length of laser beam travel; β ( r ) ( m − 1 ) is the effective volume backscatter ratio at range r , i.e., the proportion of the incident laser radiation energy that is scattered back into the solid angle subtended by the telescope receiving area at range r per unit length of laser beam travel. The quantity ( 1 / π ) Γ ( r , Ω i → Ω v ) ⋅ u L ( r ) ( s r − 1 ⋅ m − 1 ) is the volume scattering phase function, which defines the part of the radiant energy onto a unit volume of vegetative facets in the direction Ω i that is scattered in the direction Ω v within the unit solid angle. The quantity β ( r ) is the integral of the volume scattering phase function over the solid angle subtended by the telescope aperture from range r . Accordingly, we have the differential apparent reflectance ρ a p p ( r ) ( m − 1 ) from Equations (5) and (6), respectively.
ρ a p p ( r ) = J ( r ) ⋅ r 2 Φ 0 ⋅ K ( r )
ρ a p p ( r ) = P g a p ( r ) β ( r ) A T / ( π r 2 )  
Equation (9) shows how to practically calculate and interpret apparent reflectance from the received laser energy. Equation (10) shows how to theoretically model apparent reflectance from the radiative and structural characteristics of the canopy. To separate the radiative and structural information about canopy from ρ a p p ( r ) , we assume Lambertian facets and the same diffuse reflectance ρ d for all vegetative elements that contribute to the received laser energy from which ρ a p p ( r ) is calculated. Then, ( 1 / π ) Γ ( r , Ω i → Ω v ) can be approximated by G 2 ⋅ ρ d / π (see Appendix A1 for the derivation). From Equations (8) and (9), we have,
ρ a p p ( r ) = P g a p ( r ) ⋅ G 2 ⋅ u L ( r ) ⋅ ρ d
Applying the differentiation of Equation (1) to the above,
P h i t ( r ) = − ∂ P g a p ( r ) ∂ r = P g a p ( r ) ⋅ G ⋅ u L ( r )
ρ a p p ( r ) = − ∂ P g a p ( r ) ∂ r ⋅ G ⋅ ρ d
Here, P h i t ( r ) ( m − 1 ) is the laser beam interception density [53], i.e., the interception fraction by vegetative facets per unit length along the laser beam. Taking the integral over range on both sides of the Equation (13),
∫ 0 r ρ a p p ( r ′ ) d r ′ = ( 1 − P g a p ( r ) ) ⋅ G ⋅ ρ d
Thus, we can estimate P g a p over range from the integral of differential apparent reflectance over range calculated with Equation (9) from received laser energy if the G-function and leaf diffuse reflectance are known. The gap probability with range is an important function for indirect measurement of canopy structure, such as LAI, clumping index and foliage profile [17,59,60]. Calibrating lidar return intensity to apparent reflectance enables better estimates of gap probability than just counting numbers of points returned from targets along a laser beam [17].
Note that Equation (14) omits the laser beam direction and implies two assumptions: (1) the G-function is constant over range; and (2) vegetative facets are all Lambertian with the same diffuse reflectance ρ d . These two assumptions might not be true for vegetative elements traversed by a single laser beam. In practice, ρ a p p ( r ) from multiple laser shots are often averaged together (e.g., over all azimuth angles within zenith angle ranges), which will reduce the variance in estimating P g a p ( r ) for the canopy as a whole.
The foregoing discussion has not considered wavelength ( λ ) in the lidar equation and ρ a p p . From Equation (11), it is clear in the case of a bi- or multi-spectral lidar that ρ d ( λ ) of the scattering volume may be inferred from ρ a p p ( λ ) , where the bi- or multi-spectral pulse encounters the same scattering volume at the same range and from the same direction. The difference in ρ a p p ( λ ) of a scattering volume is theoretically caused only by the spectral reflectance ρ d ( λ ) , as the other terms in Equation (11) describe structural characteristics independent of wavelength. This conveniently provides a mechanism for the discrimination of different types of scattering materials based on their bi- or multi-spectral reflectance values ρ d ( λ ) .

2.1.3. Telescope Efficiency

The telescope efficiency function K ( r ) is needed by geometric laser systems using a telescope to focus the return power. It arises because the telescope is focused at infinity, and near-range objects are thus out of focus at the detector, which reduces the energy falling on the detector. The function is theoretically zero at zero range (the focal point of the telescope) and rises to unity at the range at which the focused return beam falls entirely within the detector. K ( r ) is usually omitted in the airborne lidar equation [39] because ground targets are sufficiently far enough from the instrument for K ( r ) to reach unity. For terrestrial lidar, many returns are from near-range targets, which requires including K ( r ) [30,50].

2.2. Recorded Return Waveforms and Apparent Reflectance

2.2.1. Recording Return Waveforms

The basic lidar equation above describes return laser energy ( J ) by assuming an impulse of outgoing laser energy. For a pulsing laser, the actual recording of return power over time, i.e., the return waveform, further involves: (1) the spatial distribution of outgoing laser energy (non-uniform laser beam cross-section, [61,62]); (2) the temporal shape of the outgoing laser pulse, P 0 ( t ) ( J ⋅ s − 1 ) (the spread of J 0 over a finite time); and (3) the characteristics of the detector-amplifier system of the instrument.
We shall not model the effects of the nonuniform beam cross-section on return waveforms here, because in processing terrestrial lidar data for the gap probability estimate, returns from many laser shots are usually aggregated together, which greatly reduces any variance due to the nonuniform shape of the beam cross-section.
We may express the pulse shape as a function of range, P 0 ( r ) ( J ⋅ m − 1 ), by converting the time to apparent range with r = c ( t − t p ) / 2 , where t p is the time at which the outgoing pulse peak occurs and c is the speed of light. As J 0 in Equation (4) changes to P 0 ( r ) , the term Φ 0 becomes Φ 0 ( r ) , and the return signal P ( r ) becomes the convolution of ρ a p p ( r ) and Φ 0 ( r ) :
P ( r ) = ∫ 0 r Φ 0 ( r − r ′ ) ⋅ ρ a p p   ( r ′ ) ⋅ K ( r ′ ) r ′ 2 d r ′ = Φ 0 ( r ) * ρ a p p   ( r ) ⋅ K ( r ) r 2 Φ 0 ( r ) = P 0 ( r ) A T π η s y s η a t m
As P ( r ) is altered by the detector-amplifier system, the final recorded return waveform in digital counts I ( r ) is,
I ( r ) = S R ( r ) * [ ρ a p p ( r ) ⋅ K ( r ) r 2 ]
where S R ( r ) is the result of Φ 0 ( r ) being altered by the detector-amplifier system [63].

2.2.2. Modeling Return Waveforms

To get ρ a p p ( r ) , we need to deconvolve S R ( r ) from I ( r ) . However, deconvolution is very sensitive to noise. To reduce the effects of noise on data interpretation, we modeled ρ a p p ( r ) as a sequence of Dirac delta functions marking N h ( N h ≥ 1 ) scattering clusters of vegetative facets corresponding to one or multiple return pulses in a waveform. This so-called delta-sequence model conceptually distributes vegetative facets in the j -th ( j = 1 ⋯ N h ) cluster at the range of R j inside a very thin angular voxel (i.e., Δ r → 0 ); that is, ρ a p p j ( r ) is a Dirac delta function, and a scattering cluster is a target with spatial extent of zero,
ρ a p p j ( r ) = ρ a p p j ( R j ) δ ( r − R j )
The apparent reflectance of the j -th voxel, ρ a p p j , which is the integral of differential apparent reflectance over range, is given by,
ρ a p p j = ∫ ​ ρ a p p j ( r ) d r = ∫ ​ ρ a p p j ( R j ) δ ( r − R j ) d r = ρ a p p j ( R j )
Using the delta-sequence model, I ( r ) is given by:
I ( r ) = ∑ j = 1 N h ψ j ( r ) ψ j ( r ) = S R ( r − R j ) ⋅ ρ a p p j ⋅ K ( R j ) R j 2
where ψ j ( r ) is the j-th return pulse resolvable from a return waveform. The peak intensity of ψ j ( r ) , α j is,
α j = ψ j ( R j ) = S R ( 0 ) ⋅ ρ a p p j ⋅ K ( R j ) R j 2
For a given detector-amplifier system, S R ( 0 ) is a constant, denoted as C 0 . Thus, the peak intensity of a return pulse, α j , is directly related to our calibration objective, apparent reflectance ρ a p p j of an angular voxel filled with vegetative facets as follows,
ρ a p p j = α j ⋅ R j 2 C 0 ⋅ K ( R j )
If the spatial extent of a target is not zero (violation of the Dirac delta model), ∫ ​ ρ a p p j ( r ) d r will be larger than ρ a p p j ( R j ) . In other words, the actual apparent reflectance of the angular voxel will be larger than the value calculated from the Equation (21). This underestimate of ρ a p p j due to the assumption of the Dirac delta model for clusters of vegetative facets can be compensated by correcting α j according to the shape of return pulse ψ j ( r ) and will be pursued in future.

3. Instrument and Data Preprocessing

3.1. The Dual Wavelength Echidna Lidar

The scientific objective of DWEL instrument design is to separate leaves and woody materials in forests readily in three-dimensional space using their different spectral reflectance values at near-infrared (NIR, 1064 nm) and shortwave infrared (SWIR, 1548 nm) wavelengths. At the SWIR wavelength, the laser power returned from leaves is much lower than from woody materials, such as trunks and branches, due to absorption by liquid water in leaves. In contrast, returned power from leaves and woody materials is similar at the NIR wavelength. The two infrared lasers emit unpolarized pulses with a full-width half-maximum (FWHM) of 5 ± 0.1 ns; the two laser beams are aligned coaxially to less than 1 mrad. The collimated beam diameters of the two lasers are 6 mm. The beam divergences of both lasers are 2.5 mrad for the DWEL scans presented here. The laser scanning step is 2 mrad, slightly smaller than the beam divergence, which ensures continuous coverage of the hemispheres. A third continuous-wave green marker laser is also aligned with the two infrared signal lasers; since it is readily visible, it is used to position the triple beam or mark the scan path in the laboratory. More details of the DWEL instrument design and specifications are presented in [30,38].

3.1.1. Internal Calibration Objects

Two scattering objects are fixed to the instrument to calibrate range and outgoing laser intensity. First, a fine stainless steel (removable) wire crosses the edge of the outgoing beam before it hits the scan mirror, thus scattering a small fraction of each outgoing pulse into the telescope and detectors. This allows a small part of the outgoing laser pulse to be present in the recorded signal, which assures the temporal alignment of individual waveforms and gives range precision values of one-sigma error of 4.75 cm at 1064 nm and of 2.33 cm at 1548 nm [38]. Second, a small circular Spectralon® panel with nominal reflectance of 0.99 is affixed to the case so that each mirror rotation will acquire samples of outgoing pulses from this fixed target at a fixed range. These sampled waveforms are used primarily to track drifts in laser output power that occur through the scan, but can also establish the mean outgoing pulse times of the lasers for each mirror rotation in the absence of a wire signal or refine the temporal alignment of waveforms by a wire signal.

3.1.2. Signal Recording and System Response

Figure 1 shows the mean of multiple samples of S R ( r ) given DWEL’s system response at each wavelength after background noise is removed. The pulses are normalized by peak intensity. The “ringing” response after the maximum is produced by the modulation transfer function of the combined detector-amplifier.

3.2. Preprocessing of DWEL Waveform Data

Before extraction of return pulse peaks and radiometric calibration, the raw waveforms from DWEL are preprocessed to: (1) remove background noise; (2) convert digitizer time to apparent range by aligning each waveform to the peak of the outgoing pulse using the signals from the scattering wire or internal Spectralon panel; (3) detect and correct saturated return pulses (see Section 3.2.1); (4) correct laser power drift, typically due to instrument temperature change, by scaling recorded intensities according to mean intensities observed from the internal Spectralon panel for each mirror rotation; and (5) calculate cross-covariance between waveforms and the system response function. This cross-covariance function changes the original asymmetric return pulses I ( r ) seen in each DWEL waveform to symmetric pulses. The new waveform of symmetric return pulses I p ( r ) is written as ( ⋆ denotes cross-correlation),
I p ( r ) = S R ( r ) ⋆ I ( r ) = φ ( r ) * [ ρ a p p ( r ) ⋅ K ( r ) r 2 ]
where φ ( r ) = S R ( r ) ⋆ S R ( r ) is a symmetric pulse after cross-covariance calculation. The preprocessed waveform I p ( r ) then provides the input to point cloud generation, thus avoiding pulse peak shifts due to the asymmetry of the original DWEL return pulse in later processing. An additional significant benefit of this operation is that it reduces uncorrelated noise and increases the signal-to-noise ratio prior to extraction of the signal in later processing.

Saturation Correction

Terrestrial laser scanners, in contrast to airborne scanners, will provide returns from close targets, sometimes within one meter for a placement in a forest with a dense or patchy understory, while also detecting targets at ranges of 100 m or more. This large relative variation in range provides a wide variation in return power that can exceed the limits of detector-amplifier systems (linear or nonlinear) and/or digitizers available for terrestrial scanners, and as a result, close targets can produce saturated pulse waveforms. Moreover, direct solar irradiance or specular reflectance may also produce saturated waveforms. The result may be either detector saturation, which produces an overloaded signal that persists through multiple digitizer bins or even multiple pulses, or digitizer saturation, which produces a flat-topped return pulse as the return signal exceeds the quantization range of the digitizer. In either case, the result is an unusual return pulse shape that cannot be used in calibration or to generate a scattering point with a correct apparent reflectance value.
In the DWEL instrument, detector saturation occurs in the rare case of a pulse striking an orthogonal specular target or corner reflector; normal target returns are well within the incoming power bounds of the DWEL’s detector-amplifier (Thorlabs PDA10CF) given the outgoing laser energy of DWEL. If the field of view of the telescope includes the Sun or the Sun’s aureole, the pulse may be lost completely as the detector and/or digitizer saturates or records high levels of continuous noise through the entire waveform. This situation is easy to detect, and such waveforms are identified as solar-saturated.
Digitizer saturation, however, is commonly encountered in pulses returned from near objects. Here, a “saturation correction” is employed (Figure 2). Saturation creates a flat-topped pulse as the digitizer reaches its limit; however, the side-lobe trough and secondary peak are recorded correctly. By comparing saturated and unsaturated waveforms acquired from targets with high and low reflectance at the same range, we determined empirical ratios between the magnitudes of the side lobes and the unsaturated peak. These ratios are used to generate a pulse peak that is located at the mean range of the saturated bins. This pulse is identified as a “desaturated” pulse (saturation fixed) for further processing.

4. Radiometric Calibration Procedures

4.1. Calibration Model Setup

The calibration needs to mathematically model two important components of the lidar equation, the telescope efficiency function K ( r ) and the negative exponential fall-off of return intensity with range. To set up the calibration model, we may choose between two alternatives: a physical model based on the optical design of the instrument or an empirical model designed to best fit the data. A physical model describing the returned power of the DWEL instrument was derived from first principles [30].
However, initial tests of this model with the calibration data were not satisfactory. While the general shape of the response function fits the observations, the model showed significant departures from the behavior actually observed for the instrument, especially at near range. We believe that second-order effects, such as imperfect alignment interacting with the Gaussian beam cross-section, departure of divergence from nominal specifications or unmodeled electronic effects, were responsible for this variance.
Accordingly, we used a semi-empirical model to fit the data. According to Equation (21), for each extracted point with range R and digital count intensity α , we need the constant C 0 and the function K ( r ) to calculate apparent reflectance ρ a p p . To find these unknowns, we require a collection of data points of ( R j ,   α j , ρ a p p j ) from targets of different reflectance values at multiple ranges. The quantities Rj and α j are derived from the extraction of pulse peaks; apparent reflectance ρ a p p j may be taken as the diffuse reflectance ρ d of an extended Lambertian target held perpendicular to the laser beam. The telescope efficiency function K ( r ) (see Section 2.1.3) is modeled with a generalized logistic function [64]. The empirical calibration function for the DWEL is thus:
ρ a p p ( R j ) = α j ⋅ R j b C 0 ⋅ K ( R j ) K ( r ) = 1 ( 1 + C 1 ⋅ e − C 2 ⋅ r ) C 3
where five parameters need to be estimated for the DWEL calibration, ( C 0 , C 1 , C 2 , C 3 , b ) . Of the three parameters for K ( r ) , C 1 and C 3 together affect the range at which the function approaches its asymptote of one; C 2 controls the rate at which telescope efficiency rises from zero to one in the near range.
Note that the exponent of range is now taken as a variable, b , for two reasons. First, a calibration target surface, for example a manufactured Spectralon Lambertian panel, may not provide perfectly isotropic diffuse reflectance [65]. An anisotropic target may preferentially reflect radiation into a small solid angle in the direction of the instrument’s observation (i.e., it may be specular, to some degree). Then, from Equation (2), a smaller Ω T at farther ranges would yield an integral over directions where the area scattering function, ( 1 / π ) Γ ( r , Ω i → Ω v ) , has higher values, thereby causing a larger Δ β and larger return energy J than if the target were perfectly isotropic. However, the reflectance value of a calibration panel is typically taken as a constant, for example measured by an integrating sphere, and is assumed to be isotropic in the calibration model. From Equation (5), if ρ a p p is kept constant, but J becomes larger, the exponent of range r will be smaller to compensate. Second, previous studies have also suggested that the exponent may need to accommodate electronic effects [48]. The exponent has therefore been treated as a variable in the calibration. Although the number of model parameters is large and they are not independent of each other, the calibration function can be fitted across its full range of application, thus avoiding issues of extrapolation beyond the limits of the fitting.

4.2. Calibration Data Collection

To acquire the calibration data, three panels of different reflectance values (Figure 3a) were scanned by the DWEL from a nearly perpendicular direction at 33 range locations from 0.5 m to 70 m (Table 1). The range sampling intervals were based on a provisional calibration, made at the time of commissioning, that established the general shape of the K ( r ) curve. The instrument was set in stationary mode, i.e., without the scan mirror or azimuth platform rotation (Figure 3b). The green marker laser was used to manually point the co-aligned lasers to the center of each panel at each placement (Figure 3c). The panel sizes are large enough to intercept the whole laser beam at 70 m. For each panel at each range, we collected around 150,000 waveform samples as candidates for calibration model fitting and evaluation.
The three panels included a white Spectralon panel and two foam boards painted with flat interior wall paint in light and dark gray tones derived by mixing black and white paints together. The gray panels provide unsaturated returns at near ranges. The ρ a p p value of the Spectralon panel uses the actual reflectance from the manufacturer’s specification. The ρ a p p values of the gray panels are calculated according to the Equation (7) as follows. First, we calculated the ratio of return intensities between gray panels ( J ) and the Spectralon panel ( J w ) at each range, eliminating any saturated values. Then, we averaged the product values of the ratio and Spectralon panel reflectance at all ranges to obtain the ρ a p p values of the gray panels.
Table 2 shows the ρ a p p values and the measured reflectance by a FieldSpec Pro spectrometer (Analytical Spectral Devices) of the three panels. The difference between the ρ a p p values and the measured reflectance values for the gray panels is most likely caused by two effects. First, the spectrometer measures reflectance at 0° incidence angle and 10° view angle rather than by retroreflection, and the values appeared to be underestimated slightly due to reflectance anisotropy. Moreover, the reflectance of flat wall paint may have changed with time as the paint slowly cured between the spectrometer measurements and the acquisition of calibration data (about 20 weeks).

4.3. Calibration Model Fitting

To take advantage of the difference in spectral reflectance values of leaves and those of other targets at the NIR and SWIR bands, we not only try to minimize errors in ρ a p p at individual wavelengths in calibration model fitting, but also try to ensure the two wavelengths have the same or similar relative errors in ρ a p p in order to minimize artificial variations in spectral difference due to different errors in ρ a p p at the two wavelengths. Thus, we estimate the calibration parameters of NIR and SWIR bands together in a joint calibration model. The objective error function for the fitting of this model includes both relative errors in ρ a p p from individual wavelengths and spectral constraints from both wavelengths as follows:
f ( C ˜ ) = f 1 ( C ˜ ) + f 2 ( C ˜ ) f 1 ( C ˜ ) = ∑ i = 1 N f ( ρ N i ^ − ρ N i ρ N i ) 2 + ∑ i = 1 N f ( ρ S i ^ − ρ S i ρ S i ) 2 f 2 ( C ˜ ) = v a r ( N D I ^ ) + ∑ i = 1 N f ( ρ N i ^ + ρ S i ^ − ρ N i − ρ S i ρ i + ρ S i ) 2 N D I ^ = ρ N ^ − ρ S ^ ρ N ^ + ρ S ^
where C ˜ is the vector of five calibration parameters for NIR and five parameters for SWIR; f ( C ˜ ) is the objective error function as a sum of two components, including the error f 1 ( C ˜ ) from individual wavelengths, and the spectral constraints f 2 ( C ˜ ) from both wavelengths; subscript i represents the i-th data point, and subscript N and S represent NIR and SWIR; N f is the total number of data points used in calibration fitting; ρ ^ is the apparent reflectance of panels estimated from the calibration model, while ρ is the actual apparent reflectance of panels; N D I ^ is a normalized difference index to identify the spectral difference of target reflectance between NIR and SWIR; v a r ( N D I ^ ) is the variance of NDI of data points in calibration fitting. In addition, we constrained the calibration parameters C 1 and C 3 to be equal across NIR and SWIR models. The objective error function f ( C ˜ ) has many local minima due to the high nonlinearity of the DWEL calibration model brought about by the K ( r ) function. We first used the genetic algorithm implemented in MATLAB [66] to search for initial parameter values that approach the global minimum. Then, we used the Nelder–Mead method [67] to reach the global minimum.
All of the return waveforms from the panels contain single returns at known nominal ranges. We searched maximum bins from the waveform sections around the nominal ranges and interpolated peaks using quadratic fitting of three bins around the maximum bin [68] to output intensities and ranges of these single returns. Saturated waveforms were excluded from calibration model fitting to avoid uncertainty from the saturation fix procedure. We randomly divided the remaining returns (about 24,000 samples for each range) into a training set (80 percent) and a validation set (20 percent). In the training set, return intensities were normalized by the corresponding panel reflectance to provide equivalent target reflectance values of 1.0 and then averaged together for each range to reduce noise in the data. Mean normalized intensities and mean ranges at 1064 nm and 1548 nm were paired according to panel range locations. Thirty pairs of data points from 1064 nm and 1548 nm were used to estimate the calibration parameters of 1064 nm and 1548 nm jointly by minimizing the error function f ( C ˜ ) .

5. Results and Discussion

5.1. Radiometric Calibration

5.1.1. Fitting of the Semi-Empirical Model

Table 3 provides values for the model coefficients derived from the model and procedure described in Section 4; they may be taken as examples, since recalibration will be necessary during the lifetime of the instrument.
Figure 4 and Figure 5 (first row) show the fits of the calibration functions (Equation (23)) for the two wavelengths to the training and validation data. The adjusted coefficient of determination ( R 2 ) of modeled intensity at both wavelengths for both training and validation data (Table 4) indicates that the proposed calibration function and estimated parameters predict the return intensity well. The linear regressions between measured and modeled intensity for both training and validation data (Figure 4 and Figure 5, second row) yield values of slope very close to unity, as well as very small intercepts, indicating very good fits.
The calibration function curves in Figure 4 and Figure 5 (first row), which provide the return intensity of a target with unit reflectance, increase sharply and then fall exponentially. The normalized 1064 nm return intensity peaks at ~3.5 m (Figure 4, first row) and the 1548 nm peaks at ~5 m (Figure 5, first row). The curves of K ( r ) , shown in Figure 6, rise from zero and plateau at about unity at ~10 m for the 1064 nm laser and ~15 m for 1548 nm.

5.1.2. Apparent Reflectance Error and Its Sensitivity to Intensity and Range

Because all intensities in calibration fitting and validation were normalized by corresponding panel reflectance values, the error in ρ a p p hereinafter means relative error, unless otherwise noted. The estimates of apparent reflectance ρ a p p by the calibration model from validation data show root mean squared errors (RMSE) of 0.092 at 1064 nm and 0.108 at 1548 nm (Table 4; Figure 7). The histograms of errors (Figure 7c,f) are centered near zero, which indicates little or no systematic offset in the apparent reflectance estimate.
The plots of errors in estimated ρ a p p against range for the validation dataset (Figure 7b,e) show larger and dispersed errors at very near range (<~3.5 m for 1064 nm and <~2 m for 1548 nm) and farther range (>~10 m for 1064 nm and >~20 m for 1548 nm), in contrast to smaller and less dispersed errors in between. This pattern of errors in ρ a p p over range is a combination of errors arising from range uncertainty ( Δ r ) and return intensity uncertainty ( Δ α ). We observed how Δ r and Δ α contribute to errors in ρ a p p separately with our calibration model over the range of our calibration data, 0.5 m to 70 m. We simulated δ ρ α , the relative errors in ρ a p p only due to different return intensity uncertainty levels (±15 DN) at different ranges by keeping range error at zero. Then, we simulated δ ρ r , the relative errors in ρ a p p only due to different range uncertainty levels (±15 cm) at different ranges by keeping return intensity error at zero. These relative error ranges in the simulation of Δ r and Δ α are more than three-times larger than the standard deviation of range measurements and root mean squared noise in normalized intensities.
Figure 8 presents a graphical display of the estimated relative errors calculated using the above procedure. The total relative error in ρ a p p ,   δ ρ , can be approximated by δ ρ α + δ ρ r (see Appendix A2 for this derivation). Graphs (a) and (c) in Figure 8 show the relative error in apparent reflectance ( δ ρ α ) produced by changes in return intensity ( Δ α ) of ±15 digital counts (DN) (y-axis), as the error varies with range. At very near range (<~3 m), small changes in DN produce large errors in apparent reflectance; this effect arises because the telescope efficiency K ( r ) is very low and the return signal is weak. At near range between about 2 and 10 m, the signal is much stronger, and thus, the errors produced are fairly small (green colors). Between near and far range (10 to 70 m), the exponential decrease in signal provides a smooth transition from low sensitivity of apparent reflectance with DN error to high sensitivity. At far range (70 m), the signal is sufficiently diminished by fall-off with range that errors in apparent reflectance are large given a deviation of just a few counts from true values.
In contrast, large relative error in apparent reflectance ( δ ρ r ) produced by changes in range ( Δ r ) of ±15 cm (Figure 8b,d) (y-axis) is limited to the very near range (<~5 m). The weak signal in this range provides large errors, which decrease rapidly as the telescope efficiency function K ( r ) increases the signal strength. Beyond this range, the relative error in apparent reflectance remains low.
From this analysis, we see that the error in apparent reflectance due to error in range ( δ ρ r ) dominates at near ranges, while the error due to return intensity ( δ ρ α ) dominates at far ranges. The exact range at which δ ρ α surpasses δ ρ and becomes dominant depends on the uncertainty level of return intensity given the calibration model. Thus, we see larger and more dispersed errors at very near ranges in validation data (Figure 7b,e), mainly due to range uncertainty, and at far ranges, mainly due to return intensity uncertainty. The range accuracy is therefore more critical in the near-range target calibration, while the return intensity accuracy becomes more critical in the far-range target calibration.
The problem for calibration here is that returns from far ranges have a lower signal-to-noise ratio, but their calibration is highly sensitive to return intensity uncertainty. Thus, the noise level of lidar return intensity needs to be characterized to find the range at which the reflectance uncertainty δ ρ exceeds a desirable level given the calibration model. For returns from far ranges, the apparent reflectance should be used carefully. For returns from near ranges, δ ρ could be very high if the range uncertainty is not low enough. However, lidar instruments generally give range measurements of high accuracy.

5.2. Calibration Comparison of the Two Wavelengths

The two telescope efficiency functions K ( r ) at the two wavelengths (Figure 6) suggest different optical characteristics of the two beam pathways through the DWEL instrument. As noted earlier, each pathway uses a separate wavelength-dependent divergence optic and focusing lens in the detector assembly, which can produce small differences in beam width and detector field of view. Moreover, the two laser beams may not be exactly coincident due to small errors in alignment. As a result, the two functions show different shapes.
In addition, the range exponent values b for the two wavelengths are different, but both are smaller than the theoretical value of two that applies to scattering from a perfectly diffuse surface. As noted in Section 4.1, the observed value may depart from two for a number of reasons, including slight misalignment, electronic effects in the detector-amplifier-digitizer systems and the partial specularity of target surfaces. Moreover, as a free parameter in the model inversion, the range exponent may be adjusted by the nonlinear fitting procedure to better shape the telescope efficiency function. While a more physical model grounded in instrument optics and first principles of scattering might be desirable, an empirical model will capture the data trend more accurately in the face of physical and electronic unknowns. As shown above, our calibration functions fit the data well, predicting observed intensities from calibration targets and retrieving reflectance from observed intensities with low errors.

6. Conclusions

We present a thorough derivation of the calibration objective variable, apparent reflectance, from the basic scattering lidar equation using Ross’s framework of radiation regime modeling of the vegetation canopy. As an instrument-free measure, apparent reflectance facilitates merging point data from multiple instruments, allowing assessment of the effects of angular resolution and beam divergence on structure retrieval and even providing spectral information for scanners using different laser wavelengths.
Calibration of our full-waveform, dual-wavelength, terrestrial laser scanner presents a number of challenges relevant to the next generation of terrestrial laser scanners. The need to scan from near to far range requires characterizing both telescopic effects, which reduce the near-range signal with increasing proximity due to defocusing, and saturation effects, which alter the return pulse shape of near-range scattering events. We show how to overcome these challenges by formulating a flexible calibration model, acquiring appropriate calibration data, fitting the model with a constraint providing spectral consistency and testing the results and the sensitivity of errors to uncertainties in range and intensity. We also provide solutions to the problems of saturated returns, slow change in laser output pulse energy and variance in the timing of laser pulse emissions.
In addition, dual- or multi-wavelength data must be consistent in spectral performance. By using a semi-empirical calibration model fitted to data, it is not difficult to add a constraint that optimizes spectral “flatness” with range using a Lambertian target. This step is particularly useful for the DWEL, since the laser wavelengths are chosen specifically for their ability to separate hits of water-bearing leaves from hits of the dry bark of trunks and branches and dry ground surfaces.
The RMSE values (relative errors) of apparent reflectance from our calibration procedure, 8.1 percent for 1064 nm and 6.4 percent for 1548 nm, show that the parameterized model accurately converts lidar return intensities in digital counts to apparent reflectance. This calibration model can apply to almost any terrestrial lidar instrument using a telescope to focus the return power. A sensitivity analysis shows that the apparent reflectance error from radiometric calibration is dominated by range errors at near ranges, but by return intensity errors at far ranges.
While the calibration of a terrestrial lidar is by nature more difficult than that of an airborne lidar, one advantage is that controlled laboratory or field calibration measurements can be readily designed and executed. If stationary operation can fix the beam on the target panel, it is easy to acquire pulses at measured ranges in a long corridor or outdoor environment. If stationary operation is not possible, only short scan segments crossing the target need to be acquired. Calibration is also aided by having targets of different reflectance, so that darker panels can provide unsaturated signals at ranges of the greatest returns (e.g., 10 to 12 m for DWEL). While our painted panels functioned well, a set of Lambertian panels with well-characterized diffuse reflectance properties ranging from light to dark would be desirable.
The next step is to use calibrated data to retrieve forest structural parameters with the new dual-wavelength data, following the pathways pioneered with the heritage Echidna Validation Instrument, but extending them to new information from the Dual-Wavelength Echidna lidar. These are subjects of papers now in preparation.

Acknowledgments

The authors would like to thank Jasmine Muir, Schayne Lees and Yiming Chen for their assistance in instrument operation. This research was supported by the U.S. National Science Foundation under Grant Number MRI-0923389. Echidna is a registered trademark of the Commonwealth Scientific and Industrial Research Organisation (CSIRO).

Author Contributions

Zhan Li and Alan H. Strahler designed the experiment. Zhan Li wrote the paper, with contributions from David L. B. Jupp, Alan H. Strahler and Crystal B. Schaaf. Glenn Howe, Kuravi Hewawasam, Ewan S. Douglas, Supriya Chakrabarti and Timothy A. Cook built the instrument. Zhan Li, Glenn Howe, Kuravi Hewawasam, Ian Paynter, Edward J. Saenz and Michael Schaefer carried out the experiment. Zhan Li analyzed the data. All authors reviewed and revised the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

A1. Area Scattering Phase Function for Lambertian Facets of the Same Diffuse Reflectance

Let g L ( r , Ω L ) ( s r − 1 ) be the probability function of leaf angle distribution at a location designated by range r along a laser beam of interest; the Ross G-function is [56],
G ( r , Ω ) = 1 2 π ∫ Ω 1 g L ( r , Ω L ) | cos 〈 Ω , Ω L 〉 | d Ω L
where the integral interval Ω 1 defines all of the possible leaf normal directions within an elemental volume at range r . Let ( 1 / π ) γ L ( Ω L , Ω i → Ω v ) be the leaf scattering phase function that defines the part of irradiance in the direction Ω i , which is scattered from the leaf unit area perpendicular to the direction Ω L to the unit solid angle around the direction Ω v [56]. The area scattering phase function is derived as,
1 π Γ ( r , Ω i → Ω v ) = ∫ Ω 1 g L ( r , Ω L ) 2 π γ L ( Ω L , Ω i → Ω v ) π | cos 〈 Ω i , Ω L 〉 | | cos〈 Ω v , Ω L 〉 | d Ω L
If all vegetative facets are Lambertian and have the same diffuse reflectance ρ d , since lidar instruments make observations in the backscattering direction, i.e., Ω v = Ω i − π = Ω , then we have,
1 π Γ ( r , Ω i → Ω v ) = ρ d π ∫ Ω 1 g L ( r , Ω L ) 2 π | cos〈 Ω , Ω L 〉 | 2 d Ω L ≈ ρ d π G ( r , Ω ) 2
The two G-functions in the area scattering phase function correspond to the two cosine terms in the integral. One cosine term ( | cos〈 Ω i , Ω L 〉 | ) accounts for the projected area that intercepts incoming irradiance per unit vegetative facet area. The other cosine term ( | cos 〈 Ω v , Ω L 〉 | ) accounts for the decrease of radiant intensity ( w ⋅ s r − 1 ) with the cosine, i.e., projected area for a Lambertian surface. Furthermore, we assume the leaf angle distribution is constant within the canopy and so is the G-function and area scattering function.
1 π Γ ( Ω i → Ω v ) ≈ ρ d π G ( Ω ) 2

A2. Error in Apparent Reflectance from Two Sources: Range and Return Intensity

According to Equation (23), at a given range r ,
α = C 0 K ( r ) r b ⋅ ρ a p p ρ a p p = r b C 0 K ( r ) ⋅ α
We omit the subscript j here for clarity. If the absolute errors in range r and return intensity α are Δ r and Δ α , respectively:
Δ ρ r = ( r + Δ r ) b C 0 K ( r + Δ r ) ⋅ α − ρ a p p
Δ ρ α = r b C 0 K ( r ) ⋅ ( α + Δ α ) − ρ a p p = r b C 0 K ( r ) ⋅ Δ α = ρ a p p α ⋅ Δ α
where Δ ρ r and Δ ρ α are errors in the apparent reflectance estimate due to Δ r and Δ α separately. Notice that Δ ρ r and Δ ρ α both change with range. Now, the total absolute error in apparent reflectance, Δ ρ due to Δ ρ r and Δ ρ α , together is:
Δ ρ = ( r + Δ r ) b C 0 K ( r + Δ r ) ⋅ ( α + Δ α ) − ρ a p p
Furthermore, according to Equations (A6) and (A7), we have:
( r + Δ r ) b C 0 K ( r + Δ r ) = Δ ρ r + ρ a p p α Δ α = α ⋅ Δ ρ α ρ a p p
Combining Equation (A8) and (A9),
Δ ρ = Δ ρ r + ρ a p p α ⋅ ( α + Δ ρ α ⋅ α ρ a p p ) − ρ a p p = Δ ρ r + Δ ρ α + Δ ρ r ⋅ Δ ρ α ρ a p p
The relative errors in estimated apparent reflectance due to Δ r and Δ α separately are denoted as δ ρ r = Δ ρ r / ρ a p p and δ ρ α = Δ ρ α / ρ a p p . Then, the total relative error in the apparent reflectance estimate due to Δ r and Δ α together is,
δ ρ = Δ ρ ρ a p p = 1 ρ a p p ( δ ρ r ⋅ ρ a p p + δ ρ α ⋅ ρ a p p + δ ρ r ⋅ ρ a p p ⋅ δ ρ α ⋅ ρ a p p ρ a p p ) = δ ρ r + δ ρ α + δ ρ r ⋅ δ ρ α
According to Figure 8, at near range, δ ρ r is large, while δ ρ α is very small. At far range, it is the opposite. Thus, δ ρ r ⋅ δ ρ α is small for any given range, and we can have,
δ ρ ≈ δ ρ r + δ ρ α

References

  1. Lim, K.; Treitz, P.; Wulder, M.A.; St-Onge, B.; Flood, M. LiDAR remote sensing of forest structure. Prog. Phys. Geogr. 2003, 27, 88–106. [Google Scholar] [CrossRef]
  2. Dubayah, R.O.; Drake, J.B. Lidar Remote Sensing for Forestry. J. For. 2000, 98, 44–46. [Google Scholar]
  3. Lefsky, M.A.; Cohen, W.B.; Parker, G.G.; Harding, D.J. Lidar Remote Sensing for Ecosystem Studies. Bioscience 2002, 52, 19–30. [Google Scholar] [CrossRef]
  4. Wulder, M.A.; Coops, N.C.; Hudak, A.T.; Morsdorf, F.; Nelson, R.F.; Newnham, G.J.; Vastaranta, M. Status and prospects for LiDAR remote sensing of forested ecosystems. Can. J. Remote Sens. 2013, 39, S1–S5. [Google Scholar] [CrossRef]
  5. Baldocchi, D.D.; Wilson, K.B.; Gu, L. How the environment, canopy structure and canopy physiological functioning influence carbon, water and energy fluxes of a temperate broad-leaved deciduous forest—An assessment with the biophysical model CANOAK. Tree Physiol. 2002, 22, 1065–1077. [Google Scholar] [CrossRef] [PubMed]
  6. Stark, S.C.; Leitold, V.; Wu, J.L.; Hunter, M.O.; de Castilho, C.V.; Costa, F.R.C.; McMahon, S.M.; Parker, G.G.; Shimabukuro, M.T.; Lefsky, M.A.; et al. Amazon forest carbon dynamics predicted by profiles of canopy leaf area and light environment. Ecol. Lett. 2012, 15, 1406–1414. [Google Scholar] [CrossRef] [PubMed]
  7. Strahler, A.H.; Jupp, D.L.B.; Woodcock, C.E.; Schaaf, C.L.; Yao, T.; Zhao, F.; Yang, X.; Lovell, J.L.; Culvenor, D.S.; Newnham, G.J.; et al. Retrieval of forest structural parameters using a ground-based lidar instrument. Can. J. Remote Sens. 2008, 34, S426–S440. [Google Scholar] [CrossRef]
  8. Lovell, J.L.; Jupp, D.L.B.; Newnham, G.J.; Culvenor, D.S. Measuring tree stem diameters using intensity profiles from ground-based scanning lidar from a fixed viewpoint. ISPRS J. Photogramm. Remote Sens. 2011, 66, 46–55. [Google Scholar] [CrossRef]
  9. Thies, M.; Spiecker, H. Evaluation and future prospects of terrestrial laser scanning for standardized forest inventories. Forest 2004, 2, 192–197. [Google Scholar]
  10. Huang, H.; Li, Z.; Gong, P.; Cheng, X.; Clinton, N.; Cao, C.; Ni, W.; Wang, L. Automated Methods for Measuring DBH and Tree Heights with a Commercial Scanning Lidar. Photogramm. Eng. Remote Sens. 2011, 77, 219–227. [Google Scholar] [CrossRef]
  11. Yao, T.; Yang, X.; Zhao, F.; Wang, Z.; Zhang, Q.; Jupp, D.L.B.; Lovell, J.L.; Culvenor, D.S.; Newnham, G.J.; Ni-Meister, W.; et al. Measuring forest structure and biomass in New England forest stands using Echidna ground-based lidar. Remote Sens. Environ. 2011, 115, 2965–2974. [Google Scholar] [CrossRef]
  12. Murphy, G. Determining Stand Value and Log Product Yields Using Terrestrial Lidar and Optimal Bucking: A Case Study. J. For. 2008, 106, 317–324. [Google Scholar]
  13. Raumonen, P.; Kaasalainen, M.; Åkerblom, M.; Kaasalainen, S.; Kaartinen, H.; Vastaranta, M.; Holopainen, M.; Disney, M.I.; Lewis, P. Fast Automatic Precision Tree Models from Terrestrial Laser Scanner Data. Remote Sens. 2013, 5, 491–520. [Google Scholar] [CrossRef]
  14. Béland, M.; Widlowski, J.-L.; Fournier, R.A.; Côté, J.-F.; Verstraete, M.M. Estimating leaf area distribution in savanna trees from terrestrial LiDAR measurements. Agric. For. Meteorol. 2011, 151, 1252–1266. [Google Scholar] [CrossRef]
  15. Korhonen, L.; Korpela, I.; Heiskanen, J.; Maltamo, M. Airborne discrete-return LIDAR data in the estimation of vertical canopy cover, angular canopy closure and leaf area index. Remote Sens. Environ. 2011, 115, 1065–1080. [Google Scholar] [CrossRef]
  16. Hilker, T.; Leeuwen, M.; Coops, N.C.; Wulder, M.A.; Newnham, G.J.; Jupp, D.L. B.; Culvenor, D.S. Comparing canopy metrics derived from terrestrial and airborne laser scanning in a Douglas-fir dominated forest stand. Trees 2010, 24, 819–832. [Google Scholar] [CrossRef]
  17. Jupp, D.L.B.; Culvenor, D.S.; Lovell, J.L.; Newnham, G.J.; Strahler, A.H.; Woodcock, C.E. Estimating forest LAI profiles and structural parameters using a ground-based laser called “Echidna®”. Tree Physiol. 2009, 29, 171–181. [Google Scholar] [CrossRef] [PubMed]
  18. Zhao, F.; Strahler, A.H.; Schaaf, C.L.; Yao, T.; Yang, X.; Wang, Z.; Schull, M.A.; Román, M.O.; Woodcock, C.E.; Olofsson, P.; et al. Measuring gap fraction, element clumping index and LAI in Sierra Forest stands using a full-waveform ground-based lidar. Remote Sens. Environ. 2012, 125, 73–79. [Google Scholar] [CrossRef]
  19. Dassot, M.; Constant, T.T.; Fournier, M. The use of terrestrial LiDAR technology in forest science: Application Fields, Benefits and Challenges. Ann. For. Sci. 2011, 68, 959–974. [Google Scholar] [CrossRef]
  20. Wulder, M.A.; White, J.C.; Nelson, R.F.; Næsset, E.; Ørka, H.O.; Coops, N.C.; Hilker, T.; Bater, C.W.; Gobakken, T. Lidar sampling for large-area forest characterization: A Review. Remote Sens. Environ. 2012, 121, 196–209. [Google Scholar] [CrossRef]
  21. Hancock, S.; Essery, R.; Reid, T.; Carle, J.; Baxter, R.; Rutter, N.; Huntley, B. Characterising forest gap fraction with terrestrial lidar and photography: An Examination of Relative Limitations. Agric. For. Meteorol. 2014, 189–190, 105–114. [Google Scholar] [CrossRef] [Green Version]
  22. Lefsky, M.A.; Cohen, W.B.; Acker, S.A.; Parker, G.G.; Spies, T.A.; Harding, D.J. Lidar Remote Sensing of the Canopy Structure and Biophysical Properties of Douglas-Fir Western Hemlock Forests. Remote Sens. Environ. 1999, 70, 339–361. [Google Scholar] [CrossRef]
  23. Drake, J.B.; Dubayah, R.O.; Clark, D.B.; Knox, R.G.; Blair, J.B.; Hofton, M.A.; Chazdon, R.L.; Weishampel, J.F.; Prince, S. Estimation of tropical forest structural characteristics using large-footprint lidar. Remote Sens. Environ. 2002, 79, 305–319. [Google Scholar] [CrossRef]
  24. Donoghue, D.N.M.; Watt, P.J.; Cox, N.J.; Wilson, J. Remote sensing of species mixtures in conifer plantations using LiDAR height and intensity data. Remote Sens. Environ. 2007, 110, 509–522. [Google Scholar] [CrossRef]
  25. Korpela, I.; Ørka, H.O.; Hyyppä, J.; Heikkinen, V.; Tokola, T. Range and AGC normalization in airborne discrete-return LiDAR intensity data for forest canopies. ISPRS J. Photogramm. Remote Sens. 2010, 65, 369–379. [Google Scholar] [CrossRef]
  26. Wing, B.M.; Ritchie, M.W.; Boston, K.; Cohen, W.B.; Olsen, M.J. Individual snag detection using neighborhood attribute filtered airborne lidar data. Remote Sens. Environ. 2015, 163, 165–179. [Google Scholar] [CrossRef]
  27. Gaulton, R.; Danson, F.M.; Ramirez, A.F.; Gunawan, O. The potential of dual-wavelength laser scanning for estimating vegetation moisture content. Remote Sens. Environ. 2013, 132, 32–39. [Google Scholar] [CrossRef]
  28. Eitel, J.U.H.; Magney, T.S.; Vierling, L.A.; Dittmar, G. Assessment of crop foliar nitrogen using a novel dual-wavelength laser system and implications for conducting laser-based plant physiology. ISPRS J. Photogramm. Remote Sens. 2014, 97, 229–240. [Google Scholar] [CrossRef]
  29. Kashani, A.; Olsen, M.; Parrish, C.; Wilson, N. A Review of LIDAR Radiometric Processing: From Ad Hoc Intensity Correction to Rigorous Radiometric Calibration. Sensors 2015, 15, 28099–28128. [Google Scholar] [CrossRef] [PubMed]
  30. Douglas, E.S.; Martel, J.; Li, Z.; Howe, G.A.; Hewawasam, K.; Marshall, R.A.; Schaaf, C.L.; Cook, T.A.; Newnham, G.J.; Strahler, A.H.; et al. Finding leaves in the forest: The Dual-Wavelength Echidna Lidar. Geosci. Remote Sens. Lett. IEEE 2015, 12, 776–780. [Google Scholar] [CrossRef]
  31. Danson, F.M.; Gaulton, R.; Armitage, R.P.; Disney, M.I.; Gunawan, O.; Lewis, P.; Pearson, G.; Ramirez, A.F. Developing a dual-wavelength full-waveform terrestrial laser scanner to characterize forest canopy structure. Agric. For. Meteorol. 2014, 198–199, 7–14. [Google Scholar] [CrossRef]
  32. Hakala, T.; Suomalainen, J.; Kaasalainen, S.; Chen, Y. Full waveform hyperspectral LiDAR for terrestrial laser scanning. Opt. Express 2012, 20, 7119. [Google Scholar] [CrossRef] [PubMed]
  33. Tan, S.; Narayanan, R.M. Design and Performance of a Multiwavelength Airborne Polarimetric Lidar for Vegetation Remote Sensing. Appl. Opt. 2004, 43, 2360–2368. [Google Scholar] [CrossRef] [PubMed]
  34. Gong, W.; Song, S.; Zhu, B.; Shi, S.; Li, F.; Cheng, X. Multi-wavelength canopy LiDAR for remote sensing of vegetation: Design and System Performance. ISPRS J. Photogramm. Remote Sens. 2012, 69, 1–9. [Google Scholar]
  35. Woodhouse, I.H.; Nichol, C.; Sinclair, P.; Jack, J.; Morsdorf, F.; Malthus, T.; Patenaude, G. A Multispectral Canopy LiDAR Demonstrator Project. Geosci. Remote Sens. Lett. IEEE 2011, 8, 839–843. [Google Scholar] [CrossRef]
  36. Douglas, E.S.; Strahler, A.H.; Martel, J.; Cook, T.; Mendillo, C.; Marshall, R.; Chakrabarti, S.; Schaaf, C.L.; Woodcock, C.E.; Li, Z.; et al. DWEL: A Dual-Wavelength Echidna Lidar for Ground-Based Forest Scanning. In Proceedings of the Geoscience and Remote Sensing Symposium (IGARSS), 2012 IEEE International, Munich, Germany, 26 July 2012; pp. 4998–5001.
  37. Eitel, J.U.H.; Vierling, L.A.; Long, D.S. Simultaneous measurements of plant structure and chlorophyll content in broadleaf saplings with a terrestrial laser scanner. Remote Sens. Environ. 2010, 114, 2229–2237. [Google Scholar] [CrossRef]
  38. Howe, G.A.; Hewawasam, K.; Douglas, E.S.; Martel, J.; Li, Z.; Strahler, A.H.; Schaaf, C.; Cook, T.A.; Chakrabarti, S. Capabilities and Performance of DWEL, the Dual-Wavelength Echidna Lidar. J. Appl. Remote Sens. 2015. In press. [Google Scholar]
  39. Wagner, W. Radiometric calibration of small-footprint full-waveform airborne laser scanner measurements: Basic Physical Concepts. ISPRS J. Photogramm. Remote Sens. 2010, 65, 505–513. [Google Scholar] [CrossRef]
  40. Höfle, B.; Pfeifer, N. Correction of laser scanning intensity data: Data and Model-Driven Approaches. ISPRS J. Photogramm. Remote Sens. 2007, 62, 415–433. [Google Scholar] [CrossRef]
  41. Roncat, A.; Morsdorf, F.; Briese, C.; Wagner, W.; Pfeifer, N. Laser Pulse Interaction with Forest Canopy: Geometric and Radiometric Issues. For. Appl. Airborne Laser Scanning Concepts Case Stud. 2014, 27, 19–41. [Google Scholar]
  42. Wagner, W.; Ullrich, A.; Ducic, V.; Melzer, T.; Studnicka, N. Gaussian decomposition and calibration of a novel small-footprint full-waveform digitising airborne laser scanner. ISPRS J. Photogramm. Remote Sens. 2006, 60, 100–112. [Google Scholar] [CrossRef]
  43. Kaasalainen, S.; Hyyppä, H.; Kukko, A.; Litkey, P.; Ahokas, E.; Hyyppä, J.; Lehner, H.; Jaakkola, A.; Suomalainen, J.; Akujarvi, A.; et al. Radiometric Calibration of LIDAR Intensity With Commercially Available Reference Targets. Geosci. Remote Sens. IEEE Trans. 2009, 47, 588–598. [Google Scholar] [CrossRef]
  44. Roncat, A.; Briese, C.; Jansa, J.; Pfeifer, N. Radiometrically Calibrated Features of Full-Waveform Lidar Point Clouds Based on Statistical Moments. Ieee Geosci. Remote Sens. Lett. 2014, 11, 549–553. [Google Scholar] [CrossRef]
  45. Kaasalainen, S.; Pyysalo, U.; Krooks, A.; Vain, A.; Kukko, A.; Hyyppä, J.; Kaasalainen, M. Absolute Radiometric Calibration of ALS Intensity Data: Effects on Accuracy and Target Classification. Sensors 2011, 11, 10586–10602. [Google Scholar] [CrossRef] [PubMed]
  46. Briese, C.; Pfennigbauer, M.; Lehner, H.; Ullrich, A.; Wagner, W.; Pfeifer, N. Radiometric calibration of multi-wavelength airborne laser scanning data. ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci. 2012, I-7, 335–340. [Google Scholar] [CrossRef]
  47. Béland, M.; Baldocchi, D.D.; Widlowski, J.-L.; Fournier, R.A.; Verstraete, M.M. On seeing the wood from the leaves and the role of voxel size in determining leaf area distribution of forests with terrestrial LiDAR. Agric. For. Meteorol. 2014, 184, 82–97. [Google Scholar] [CrossRef]
  48. Yang, X.; Strahler, A.H.; Schaaf, C.L.; Jupp, D.L.B.; Yao, T.; Zhao, F.; Wang, Z.; Culvenor, D.S.; Newnham, G.J.; Lovell, J.L.; et al. Three-dimensional forest reconstruction and structural parameter retrievals using a terrestrial full-waveform lidar instrument (Echidna®). Remote Sens. Environ. 2013, 135, 36–51. [Google Scholar] [CrossRef]
  49. Hartzell, P.; Glennie, C.; Finnegan, D.C. Empirical Waveform Decomposition and Radiometric Calibration of a Terrestrial Full-Waveform Laser Scanner. Geosci. Remote Sens. IEEE Trans. 2015, 53, 162–172. [Google Scholar] [CrossRef]
  50. Pfennigbauer, M.; Ullrich, A. Improving quality of laser scanning data acquisition through calibrated amplitude and pulse deviation measurement. Proc. SPIE 2010, 7684. [Google Scholar] [CrossRef]
  51. Pfeifer, N.; Höfle, B.; Briese, C.; Rutzinger, M.; Haring, A. Analysis of the backscattered energy in terrestrial laser scanning data. In Proceedings of the XXIth ISPRS Congress, Silk Road for Information from Imagery, Beijing, China, 8 July 2008; Volume 37, pp. 1045–1052.
  52. Kaasalainen, S.; Krooks, A.; Kukko, A.; Kaartinen, H. Radiometric Calibration of Terrestrial Laser Scanners with External Reference Targets. Remote Sens. 2009, 1, 144–158. [Google Scholar] [CrossRef]
  53. Ni-Meister, W.; Jupp, D.L. B.; Dubayah, R.O. Modeling lidar waveforms in heterogeneous and discrete canopies. Geosci. Remote Sensing IEEE Trans. 2001, 39, 1943–1958. [Google Scholar] [CrossRef]
  54. Sun, G.; Ranson, K.J. Modeling lidar returns from forest canopies. Geosci. Remote Sensing IEEE Trans. 2000, 38, 2617–2626. [Google Scholar]
  55. Measures, R.M. Laser-remote-sensor equations. In Laser Remote Sensing: Fundamentals and Applications; Wiley-Interscience: New York, NY, USA, 1984; pp. 237–243. [Google Scholar]
  56. Ross, J. The Radiation Regime and Architecture of Plant Stands; Lieth, H., Ed.; Dr W. Junk Publishers: The Hague, The Netherlands, 1981. [Google Scholar]
  57. Hancock, S. Understanding the Measurement of Forests with Waveform Lidar. Ph.D. Thesis, University College London, London, UK, 2010. [Google Scholar]
  58. Parkin, D.A.; Jupp, D.L.B.; Poropat, G.V; Lovell, J.L. Lidar System and Method. CSIRO PCT TW6959. 9 February 2001. [Google Scholar]
  59. Welles, J.M.; Cohen, S. Canopy structure measurement by gap fraction analysis using commercial instrumentation. J. Exp. Bot. 1996, 47, 1335–1342. [Google Scholar] [CrossRef]
  60. Chen, J.M.; Cihlar, J. Plant canopy gap-size analysis theory for improving optical measurements of leaf-area index. Appl. Opt. 1995, 34, 6211–6222. [Google Scholar] [CrossRef] [PubMed]
  61. Harding, D.J.; Lefsky, M.A.; Parker, G.G.; Blair, J.B. Laser altimeter canopy height profiles: Methods and Validation for Closed-Canopy, Broadleaf Forests. Remote Sens. Environ. 2001, 76, 283–297. [Google Scholar] [CrossRef]
  62. Jutzi, B.; Eberle, B.; Stilla, U. Estimation and measurement of backscattered signals from pulsed laser radar. Proc. SPIE 2003, 4885. [Google Scholar] [CrossRef]
  63. Roncata, A.; Lehnera, H.; Briesea, C. Laser pulse variations and their influence on radiometric calibration of full-waveform laser scanner data. ISPRS-International Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2011, 3812, 37–42. [Google Scholar] [CrossRef]
  64. Richards, F.J. A Flexible Growth Function for Empirical Use. J. Exp. Bot. 1959, 10, 290–301. [Google Scholar] [CrossRef]
  65. Bhandari, A.; Hamre, B.; Frette, Ø.; Zhao, L.; Stamnes, J.J.; Kildemo, M. Bidirectional reflectance distribution function of Spectralon white reflectance standard illuminated by incoherent unpolarized and plane-polarized light. Appl. Opt. 2011, 50, 2431–2442. [Google Scholar] [CrossRef] [PubMed]
  66. MATLAB Genetic Algorithm. Avaliable online: http://www.mathworks.com/discovery/genetic-algorithm.html (accessed on 20 January 2015).
  67. MATLAB Fminsearch Algorithm. Avaliable online: http://www.mathworks.com/help/matlab/math/optimizing-nonlinear-functions.html#bsgpq6p-11 (accessed on 20 January 2015).
  68. Hancock, S.; Armston, J.; Li, Z.; Gaulton, R.; Lewis, P.; Disney, M.I.; Danson, F.M.; Strahler, A.H.; Schaaf, C.L.; Anderson, K.; et al. Waveform lidar over vegetation: An Evaluation of Inversion Methods for Estimating Return Energy. Remote Sens. Environ. 2014, 164, 208–224. [Google Scholar] [CrossRef] [Green Version]
Sensors 16 00313 g001 1024
Figure 1. DWEL system response pulse. The two pulse peaks at the two wavelengths are aligned.
Figure 1. DWEL system response pulse. The two pulse peaks at the two wavelengths are aligned.
Sensors 16 00313 g001
Sensors 16 00313 g002 1024
Figure 2. Example of saturated pulse and saturation correction.
Figure 2. Example of saturated pulse and saturation correction.
Sensors 16 00313 g002
Sensors 16 00313 g003 1024
Figure 3. Calibration data collection. (a) Three panels for calibration, from top to bottom: painted light gray panel, white Spectralon panel and painted dark gray panel; (b) DWEL was set up in stationary mode with the laser pointing along the measuring tape laid out on the floor; panels were placed at 33 range locations along the tape; (c) the green laser was used to point the infrared lasers to the center of the panels.
Figure 3. Calibration data collection. (a) Three panels for calibration, from top to bottom: painted light gray panel, white Spectralon panel and painted dark gray panel; (b) DWEL was set up in stationary mode with the laser pointing along the measuring tape laid out on the floor; panels were placed at 33 range locations along the tape; (c) the green laser was used to point the infrared lasers to the center of the panels.
Sensors 16 00313 g003
Sensors 16 00313 g004 1024
Figure 4. Estimation and validation of calibration of 1064 nm data. In both rows, the left column (a,c) shows the calibration function as fitted to training data, and the right column (b,d) shows the fit to the validation data. First row (a,b): measured and modeled intensity normalized by reflectance. Second row (c,d): scatter plots of measured against modeled intensity. The vertical error bars in (a) and horizontal error bars in (c) are one standard deviation of measured intensities normalized by reflectance.
Figure 4. Estimation and validation of calibration of 1064 nm data. In both rows, the left column (a,c) shows the calibration function as fitted to training data, and the right column (b,d) shows the fit to the validation data. First row (a,b): measured and modeled intensity normalized by reflectance. Second row (c,d): scatter plots of measured against modeled intensity. The vertical error bars in (a) and horizontal error bars in (c) are one standard deviation of measured intensities normalized by reflectance.
Sensors 16 00313 g004
Sensors 16 00313 g005 1024
Figure 5. Estimation and validation of calibration of 1548 nm data. In both rows, the left column (a,c) shows the calibration function as fitted to training data, and the right column (b,d) shows the fit to the validation data. First row (a,b): measured and modeled intensity normalized by reflectance. Second row (c,d): scatter plots of measured against modeled intensity. The vertical error bars in (a) and horizontal error bars in (c) are one standard deviation of measured intensities normalized by reflectance.
Figure 5. Estimation and validation of calibration of 1548 nm data. In both rows, the left column (a,c) shows the calibration function as fitted to training data, and the right column (b,d) shows the fit to the validation data. First row (a,b): measured and modeled intensity normalized by reflectance. Second row (c,d): scatter plots of measured against modeled intensity. The vertical error bars in (a) and horizontal error bars in (c) are one standard deviation of measured intensities normalized by reflectance.
Sensors 16 00313 g005
Sensors 16 00313 g006 1024
Figure 6. Telescope efficiency K ( r ) of the two wavelengths.
Figure 6. Telescope efficiency K ( r ) of the two wavelengths.
Sensors 16 00313 g006
Sensors 16 00313 g007 1024
Figure 7. Errors in apparent reflectance. The first row (a–c) shows 1064 nm, and the second row (d–f) shows 1548 nm. The left column (a,d) is the deviation from calibration fitting with range. The middle (b,e) column is the deviation of validation points with range. The right column (c,f) is the histogram of deviations.
Figure 7. Errors in apparent reflectance. The first row (a–c) shows 1064 nm, and the second row (d–f) shows 1548 nm. The left column (a,d) is the deviation from calibration fitting with range. The middle (b,e) column is the deviation of validation points with range. The right column (c,f) is the histogram of deviations.
Sensors 16 00313 g007
Sensors 16 00313 g008 1024
Figure 8. Sensitivity of the ρ a p p estimate on errors in return intensity and range. The image color shows relative error in the ρ a p p estimate (estimate − measurement). The color map scale is unified for all images for comparison purposes, but the actual error ranges of the four images are different and given here: (a) δ ρ α at 1064 nm [−0.928, 0.928]; (b) δ ρ r at 1064 nm, [−0.226, 0.290]; (c) δ ρ α at 1548 nm, [−0.574, 0.574]; (d) δ ρ r at 1548, [−0.133, 0.154].
Figure 8. Sensitivity of the ρ a p p estimate on errors in return intensity and range. The image color shows relative error in the ρ a p p estimate (estimate − measurement). The color map scale is unified for all images for comparison purposes, but the actual error ranges of the four images are different and given here: (a) δ ρ α at 1064 nm [−0.928, 0.928]; (b) δ ρ r at 1064 nm, [−0.226, 0.290]; (c) δ ρ α at 1548 nm, [−0.574, 0.574]; (d) δ ρ r at 1548, [−0.133, 0.154].
Sensors 16 00313 g008
Table
Table 1. Range sample design.
Table 1. Range sample design.
Range (m)Range Interval (m)Measurement Positions
[0.5, 10] *0.520
(10, 15]15
(15, 40]55
(40, 70]103
* “[” or “]” means inclusive of range boundary while “(” means exclusive of range boundary.
Table
Table 2. Reflectance values of panels used in calibration.
Table 2. Reflectance values of panels used in calibration.
TargetNIR ReflectanceSWIR ReflectanceDimension (cm by cm)
Measured 1 ρ a p p   2Measured 1 ρ a p p   2
White Spectralon panel 30.990.9830.5 × 30.5
Gray Painted Panel 10.4360.5740.3490.44738.0 × 30.5
Gray Painted Panel 20.3200.4310.2450.32938.0 × 30.5
1 From the spectrometer with the illuminated probe; 2 calculated according to Equation (7); 3 manufacturer’s calibrated value.
Table
Table 3. An example set of DWEL calibration parameters.
Table 3. An example set of DWEL calibration parameters.
ParameterWavelength
1064 nm1548 nm
C05788.26581822,054.218342
C10.0003190.000319
C20.8088800.540762
C325,176.83503225,176.835032
b1.3842971.585985
Table
Table 4. Assessment of calibration fitting and validation.
Table 4. Assessment of calibration fitting and validation.
Wavelength1064 nm1548 nm
Measured vs. Modeled Intensity, Adjusted R 2 Training 10.9540.983
Validation 20.9480.964
RMSE of Apparent ReflectanceTraining 10.1080.092
Validation 20.0810.064
1 Thirty data points used for training by averaging unsaturated waveforms from three panels at 30 ranges (19,200 waveforms at each range); no good waveforms from three ranges (0.5, 1.0 and 70 m) out of 33 measured locations; 2 about 4800 waveforms used for validation. See Section 4.2 and Section 4.3 for more details.

Share and Cite

MDPI and ACS Style

Li, Z.; Jupp, D.L.B.; Strahler, A.H.; Schaaf, C.B.; Howe, G.; Hewawasam, K.; Douglas, E.S.; Chakrabarti, S.; Cook, T.A.; Paynter, I.; et al. Radiometric Calibration of a Dual-Wavelength, Full-Waveform Terrestrial Lidar. Sensors 2016, 16, 313. https://doi.org/10.3390/s16030313

AMA Style

Li Z, Jupp DLB, Strahler AH, Schaaf CB, Howe G, Hewawasam K, Douglas ES, Chakrabarti S, Cook TA, Paynter I, et al. Radiometric Calibration of a Dual-Wavelength, Full-Waveform Terrestrial Lidar. Sensors. 2016; 16(3):313. https://doi.org/10.3390/s16030313

Chicago/Turabian Style

Li, Zhan, David L. B. Jupp, Alan H. Strahler, Crystal B. Schaaf, Glenn Howe, Kuravi Hewawasam, Ewan S. Douglas, Supriya Chakrabarti, Timothy A. Cook, Ian Paynter, and et al. 2016. "Radiometric Calibration of a Dual-Wavelength, Full-Waveform Terrestrial Lidar" Sensors 16, no. 3: 313. https://doi.org/10.3390/s16030313

APA Style

Li, Z., Jupp, D. L. B., Strahler, A. H., Schaaf, C. B., Howe, G., Hewawasam, K., Douglas, E. S., Chakrabarti, S., Cook, T. A., Paynter, I., Saenz, E. J., & Schaefer, M. (2016). Radiometric Calibration of a Dual-Wavelength, Full-Waveform Terrestrial Lidar. Sensors, 16(3), 313. https://doi.org/10.3390/s16030313

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop