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

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.


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.

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 gap prq : where G is the Ross G-function, which describes the projection of a unit vegetative area in a given direction; u L prq (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 gap prq 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: In these equations, J 0 P gap prq (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 p1{πq Γ pr, Ω i Ñ Ω v q (sr´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 (sr) is the solid angle subtended by the area of the telescope aperture (A T ) of the lidar instrument from range r; K prq (dimensionless) is telescope efficiency at range r (see Section 2.1.3); and η sys and η atm (dimensionless) are transmission factors that account for energy loss due to the sensor system and atmosphere, respectively.

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 gap prq " 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 p1{πq Γ pΩ i Ñ Ω v q " ρ d {π; and (3) all of the vegetative facets together fill the whole laser beam, that is r`∆r r u L`r 1˘d r 1 " 1. The return energy from such a voxel is then: 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ȋ n Equation (3) is the proportion of the total hemispherical backscattering that is intercepted by the telescope aperture given this simplest case. Assuming η sys and η atm 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., 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, ρ app (dimensionless), to represent the overall structural and radiative characteristics of all of the vegetative facets as a whole illuminated by a laser beam. The ρ app for the general case is calculated in the same way as ρ d for the simplest case, i.e., ∆J " Φ 0¨ρ app¨K prq r 2 (5) Comparing Equation (5) to Equation (2), we also have, 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, ρ app can be modeled to derive canopy radiative and structural information. According to Equation (6), ρ app includes P gap prq, 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, ρ app 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): Thus, ρ app 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].

Physical Interpretation of Apparent Reflectance
When the thickness of the angular voxel ∆r Ñ 0 , we have the differential form of Equation (2) as: where J prq (J¨m´1) is the received laser energy from range r per unit length of laser beam travel; β prq (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 p1{πq Γ pr, Ω i Ñ Ω v q¨u L prq (sr´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 β prq 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 ρ app prq (m´1) from Equations (5) and (6), respectively.
ρ app prq " P gap prq β prq A T {`πr 2˘( 10) 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 ρ app prq, we assume Lambertian facets and the same diffuse reflectance ρ d for all vegetative elements that contribute to the received laser energy from which ρ app prq is calculated. Then, p1{πq Γ pr, Ω i Ñ Ω v q can be approximated by G 2¨ρ d {π (see Appendix A1 for the derivation). From Equations (8) and (9), we have, Applying the differentiation of Equation (1) to the above, P hit prq "´B P gap prq Br " P gap prq¨G¨u L prq (12) ρ app prq "´B P gap prq Here, P hit prq (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), Thus, we can estimate P gap 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, ρ app prq 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 gap prq for the canopy as a whole.
The foregoing discussion has not considered wavelength (λ) in the lidar equation and ρ app . From Equation (11), it is clear in the case of a bi-or multi-spectral lidar that ρ d pλq of the scattering volume may be inferred from ρ app pλq, where the bi-or multi-spectral pulse encounters the same scattering volume at the same range and from the same direction. The difference in ρ app pλq of a scattering volume is theoretically caused only by the spectral reflectance ρ d pλq, 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 pλq.

Telescope Efficiency
The telescope efficiency function K prq 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 prq is usually omitted in the airborne lidar equation [39] because ground targets are sufficiently far enough from the instrument for K prq to reach unity. For terrestrial lidar, many returns are from near-range targets, which requires including K prq [30,50].

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 ptq (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 prq (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 prq, the term Φ 0 becomes Φ 0 prq, and the return signal P prq becomes the convolution of ρ app prq and Φ 0 prq: As P prq is altered by the detector-amplifier system, the final recorded return waveform in digital counts I prq is, where S R prq is the result of Φ 0 prq being altered by the detector-amplifier system [63].

Modeling Return Waveforms
To get ρ app prq, we need to deconvolve S R prq from I prq. However, deconvolution is very sensitive to noise. To reduce the effects of noise on data interpretation, we modeled ρ app prq 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, ρ app j prq is a Dirac delta function, and a scattering cluster is a target with spatial extent of zero, The apparent reflectance of the j-th voxel, ρ app j , which is the integral of differential apparent reflectance over range, is given by, Using the delta-sequence model, I prq is given by: where ψ j prq is the j-th return pulse resolvable from a return waveform. The peak intensity of For a given detector-amplifier system, S R p0q 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 ρ app j of an angular voxel filled with vegetative facets as follows, If the spatial extent of a target is not zero (violation of the Dirac delta model), ρ app j prq dr will be larger than ρ app 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 ρ app 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 prq and will be pursued in future.

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].

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. Figure 1 shows the mean of multiple samples of S R prq 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.

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 ( ) seen in each DWEL waveform to symmetric pulses. The new waveform of symmetric return pulses ( ) is written as (⋆ denotes cross-correlation), where ( ) = ( ) ⋆ ( ) is a symmetric pulse after cross-covariance calculation. The preprocessed waveform ( ) 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-tonoise 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

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); (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 prq seen in each DWEL waveform to symmetric pulses. The new waveform of symmetric return pulses I p prq is written as (‹ denotes cross-correlation), where ϕ prq " S R prq ‹ S R prq is a symmetric pulse after cross-covariance calculation. The preprocessed waveform I p prq 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. 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.

Calibration Model Setup
The calibration needs to mathematically model two important components of the lidar equation, the telescope efficiency function ( ) 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 and digital count intensity , we need the constant 0 and the function ( ) to calculate apparent reflectance . To find these unknowns, we require a collection of data points of ( , , ) from targets of different reflectance values at multiple ranges. The quantities Rj and are derived from the extraction of pulse peaks; apparent reflectance may be taken as the diffuse reflectance of an extended Lambertian target held perpendicular

Calibration Model Setup
The calibration needs to mathematically model two important components of the lidar equation, the telescope efficiency function K prq 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 prq to calculate apparent reflectance ρ app . To find these unknowns, we require a collection of data points of´R j , α j , ρ app j¯f rom targets of different reflectance values at multiple ranges. The quantities R j and α j are derived from the extraction of pulse peaks; apparent reflectance ρ app 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 prq (see Section 2.1.3) is modeled with a generalized logistic function [64]. The empirical calibration function for the DWEL is thus: where five parameters need to be estimated for the DWEL calibration, pC 0 , C 1 , C 2 , C 3 , bq. Of the three parameters for K prq, 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, p1{πq Γ pr, Ω i Ñ Ω v q, 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 ρ app 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.

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 prq 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 ρ app value of the Spectralon panel uses the actual reflectance from the manufacturer's specification. The ρ app 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 ρ app values of the gray panels.  Table 2 shows the ρ app values and the measured reflectance by a FieldSpec Pro spectrometer (Analytical Spectral Devices) of the three panels. The difference between the ρ app 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). intensities between gray panels ( ) and the Spectralon panel ( ) 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 values of the gray panels. Table 2 shows the values and the measured reflectance by a FieldSpec Pro spectrometer (Analytical Spectral Devices) of the three panels. The difference between the 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).

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 at individual wavelengths in calibration model fitting, but also try to ensure the two wavelengths have the same or similar relative errors in in order to minimize artificial variations in spectral difference due to different errors in at the two wavelengths. Thus, we estimate the calibration parameters of

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 ρ app at individual wavelengths in calibration model fitting, but also try to ensure the two wavelengths have the same or similar relative errors in ρ app in order to minimize artificial variations in spectral difference due to different errors in ρ app 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 ρ app from individual wavelengths and spectral constraints from both wavelengths as follows: 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˘f rom individual wavelengths, and the spectral constraints f 2`C˘f rom 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 DI is a normalized difference index to identify the spectral difference of target reflectance between NIR and SWIR; var`N DI˘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 prq 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˘. 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.  (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) (23)) for the two wavelengths to the training and validation data. The adjusted coefficient of determination ( 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 (Figures 4 and 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 Figures 4 and 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 ( ), 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.  The calibration function curves in Figures 4 and 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 prq, 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. The plots of errors in estimated 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 over range is a combination of errors arising from range uncertainty ( Δ ) and return intensity uncertainty ( Δ ). We observed how Δ and Δ contribute to errors in separately with our calibration model over the range of our calibration data, 0.5 m to 70 m. We simulated , the relative errors in only due to different return intensity uncertainty levels (±15 DN) at different ranges by keeping range error at zero. Then, we simulated , the relative errors in 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 Δ 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 , , can be approximated by + (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 ( ) 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.

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 ρ app hereinafter means relative error, unless otherwise noted. The estimates of apparent reflectance ρ app 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 ρ app 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 ρ app 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 ρ app separately with our calibration model over the range of our calibration data, 0.5 m to 70 m. We simulated δρ α , the relative errors in ρ app 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 ρ app 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 ρ app , δρ, 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 prq 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 ( ) produced by changes in range (Δ ) 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 ( ) 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 ( ) 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. In contrast, large relative error in apparent reflectance (δρ r ) produced by changes in range (∆ rq 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 prq 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 q 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.

Calibration Comparison of the Two Wavelengths
The two telescope efficiency functions K prq 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.

Calibration Comparison of the Two Wavelengths
The two telescope efficiency functions ( ) 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 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. 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.

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.