CORRECTING ATTENUATION EFFECTS CAUSED BY INTERACTIONS IN THE FOREST CANOPY IN FULL-WAVEFORM AIRBORNE LASER SCANNER DATA

Full-waveform airborne laser scanning offers a great potential for v arious forestry applications. Especially applications requiring information on the vertical structure of the lower canopy parts benefit fr om the great amount of information contained in waveform data. To enable the derivation of vertical forest canopy structure, the dev lopment of suitable voxel based data analysis methods is straightforward. Beyond extracting additional 3D points, it is very promis ing to derive the voxel attributes from the digitized waveform directly. For this purpose, the differential backscatter cross sections h ave to be projected into a Cartesian voxel structure. Thereby the voxel entries represent amplitudes of the cross section and can be inter preted as a local measure for the amount of pulse reflecting matter. However, the ’history’ of each laser echo pulse is characterized by atte nu tion effects caused by reflections in higher regions of the crown. As a result, the received waveform signals within the canopy hav e a lower amplitude than it would be observed for an identical structure without the previous canopy structure interactions (Romanczy k et al., 2012). If the biophysical structure is determined from the raw waveform data, material in the lower parts of the canopy is thus und er-represented. To achieve a radiometrically correct voxel space representation the los s of signal strength caused by partial reflections on the path of a laser pulse through the canopy has to be compensated. In this paper, we present an integral approach correcting the waveform at each recorded sample. The basic idea of the procedure is to enhance the wav eform intensity values in lower parts of the canopy for portions of the pulse intensity, which have been reflected (and thus blocked) in high er parts of the canopy. The paper will discuss the developed correction method and show results from a validation both with synthetic and r eal world data.


INTRODUCTION
The characterization of the vertical forest structure is highly relevant for ecological research and for better understanding forest ecosystems.Small-footprint full-waveform laser scanning has been increasingly used in forestry applications for some years now, because it delivers measurements on physical attributes of vegetation canopy structure.Voxel based data-analysis methods have a great potential to derive biophysical forest parameters.Lindberg et al. (2012) point out that using ALS waveform data to describe the volumetric aspects of vertical vegetation structure is very promising.So far, new approaches are mostly based on the generation of 3D raster domains derived from extracted discrete 3D points (e.g.Reitberger et al. (2009), Leiterer et al. (2012), Hosoi et al. (2013)).Deriving voxel attributes from the digitized waveform directly instead of using a densified point cloud was proposed by Persson et al. (2005) and Buddenbaum et al. (2013).
We expect that these volumetric representations will facilitate the determination of parameters which could so far hardly be extracted from airborne laser scanner data, for instance vertical and horizontal forest structure, open stem area, lower tree crown delineation and classification as well as the density of understory vegetation.Moreover, a significant improvement of the precision of biomass and biodiversity parameters can be expected.Generating volumetric forest stand representations of high geometric and radiometric quality from full-waveform airborne laser scanner data requires a series of geometric and radiometric transformations of the data into a voxel structure as well as an improved understanding of the interaction between laser pulse and pulse reflecting matter.
An important factor herein are attenuation effects caused by reflections in higher regions of the crown.So far, these effects are either not considered or only considered in a very simple way.Lindberg et al. (2012) use an iterative normalization algorithm based on the Beer-Lambert law to compensate the attenuation of the signal energy as the signal propagates through a canopy.Allouis et al. (2011) correct the aggregated waveform with a logarithmic function.Furthermore, several simulation studies contribute to an improved understanding of signal-forest interactions (Sun and Ranson (2000), Liu et al. (2010), Pang et al. (2011), Romanczyk et al. (2012), Cawse-Nicholson et al. (2013)).
In order to derive the biophysical structure of a forest stand from full-waveform airborne laser scanner data, advanced methods for attenuation correction are required.Finding a suitable correction method will contribute to improve the accuracy of volumetric forest reconstructions and enable forestry applications which require accurate information in the understory of a forest canopy (Cawse-Nicholson et al., 2013).In this paper, we present an attenuation correction method based on a waveform history analysis.
The paper is divided into 5 parts: Part 2 treats the theoretical background, and part 3 describes the developed attenuation correction model.The results of the validation with both synthetic and real world data are discussed in part 4. Finally a short conclusion is given in part 5.

PHYSICAL PRINCIPLES AND RELATED WORK
The physical principle of full-waveform airborne laser scanning comprises the emission of laser pulses, the interaction of the emitted pulse with the target, and the recording of the backscattered signal, called received waveform.On one hand, the shape of the received waveform depends on several sensor parameters (shape of the laser pulse, receiver impulse function, pulse spreading), on the other hand on the backscattering characteristics of the target (Wagner et al., 2007).For spatially distributed targets, the return signal is mainly a superposition of echoes from scatterers at different ranges.If the distance between the scatterers is larger than the range resolution of the ALS system, they produce distinct echoes.Normally, the leaves and branches of vegetation are clustered at smaller distances.If the emitted laser pulse interacts with such a cluster, the original pulse is widened.All target properties, e.g.reflectivity, target size, and directionality of scattering, are summarized into one parameter: the differential backscatter cross section σ(t).The formation of the received waveform Pr(t) can be described mathematically as the result of a convolution of the emitted laser pulse Pe(t) and the differential backscatter cross section.As shown in equation 1, the received signal depends on the aperture diameter of the receiver optic Dr, the range R, and the transmitter beam width β 2 t (Wagner et al., 2006).
Depending on the chosen method, the cross section is treated as a sum of discrete values or as a continuous variable (Wagner et al., 2006).In the latter case, the differential backscatter cross section can be projected into a Cartesian voxel structure (Figure 1).The voxel entries represent the amplitudes of the corresponding differential backscatter cross section.The procedure of filling the voxel space has both geometric and radiometric aspects.The geometric aspects include an intersection of the diverging laser cone with the voxel structure and an interpolation procedure to obtain the actual voxel values.The differential backscatter cross section resulting from decovolution is in first instance only an "apparent" cross section, because the measuring process is influenced by different attenuation effects.Therefore, the amplitudes of the backscatter cross section have to undergo some radiometric correction before entering them into the voxel structure.There are two groups of attenuation effects which influence the emitted laser pulse and therefore the shape of received waveform and derived differential backscatter cross section -on the one hand the attenuation due to atmospheric effects, on the other hand the attenuation of the signal during its propagation through the canopy.
The former is treated in several studies (Höfle and Pfeifer (2007), Wagner ( 2010)) and does not require a further discussion here.The latter issue is analyzed in more detail in the following section.After the radiometric correction, the voxel entries can be interpreted as a local measure for the amount of pulse reflecting matter.
In their simulation study, Cawse-  derive range, amplitude and standard deviation for each peak, and developed a correction factor based on the results of the Gaussian decomposition.Our correction method, that we applied to both simulated and real data in this study, is similar to their approach, even though it is not based on the detailed consideration of single photons.In contrast to their approach treating the attenuation theoretically in a simulation, our investigation starts from a practical point of view (unknown tree geometry and reflectance).Our objective is to develop a correction method, which can be easily applied to real full-waveform forestry data sets.Thereby, the correction term has to be derived from the digitized pulse echo itself.Furthermore, our attenuation correction is derived from the complete recorded signal, not applying decomposition methods.

METHOD
During the propagation of the emitted laser pulse through a forest canopy, it interacts with different components (leaves, branches, forest floor).Depending on the material, the number of incoming photons is reduced due to reflection, transmission, absorption and scattering processes.Therefore, less photons are available for subsequent interactions in lower parts of the canopy.As a result, the received waveform signal within the canopy has a lower amplitude than it would be observed for an identical structure without the previous (upper) canopy structure interactions (Romanczyk et al., 2013).If the biophysical structure is determined from the waveform raw data, the material in the lower parts of the canopy is thus not estimated correctly.The following considerations require that the attenuation due to atmospheric effects has already been corrected.

Attenuation Model
First of all we need to describe the attenuation during the propagation of the laser pulse through the canopy with a suitable model.This can be done either segment-wise in discrete steps or integral for each sample.In (Richter et al., 2014), we presented an attenuation correction method based on a discrete attenuation model.The model assumes that the pulse intensity is decreased at each cluster of scatterers by a certain factor which depends on the portion of the laser pulse reflected at the previous interactions.In the approach presented in this paper, we assume that the pulse intensity is decreased continuously during the propagation, resulting in the extension of the discrete model to an integral one.Please note that this approach modifies the shape of the backscattered signal and influences the geometric component of the signal.Strictly speaking, this is only valid for volumetric scatterers like leaves and branches, not for extended targets like the forest floor.In the following, the difference of the two attenuation models is first illustrated by means of our synthetic data set.
The correction methods based on both models are discussed in 4. taking into account synthtic data as well as real world data.

Synthetic Data Set
For validation purposes, we created a synthetic example (Figure 2) consisting of two cases: the undisturbed propagation of the emitted laser pulse to the ground without any interactions (left) which serves as reference, and the propagation through a canopy with interactions at three clusters and the forest floor (right).We designed a synthetic emitted waveform and two synthetic differential backscatter cross sections, where the spatial expansion of the tree layers is considered.Please note that the cross sections are true differential backscatter cross sections directly correlated to the amount of pulse reflecting matter.An ideal received waveform without attenuation effects is obtained by a convolution of emitted waveform and cross section.In the next step, the attenuation is taken into account by decreasing the received waveform with a linear attenuation factor ai.The factor depends on the portion of the laser pulse which is reflected at each interaction (discrete model) or waveform digitization sample (integral model).
The reflected portion can be estimated from the relation between reference cross section and canopy cross section.For the discrete attenuation model (Figure 3), the attenuation factor is calculated from the ratio between the area under the reference cross section (Figure 2, blue) and the area under the cross section for each interaction A1 − Ai (Figure 2, green).
For the extension to the integral model (Figure 4), the areas A1 − Ai are superseded by the amplitudes csi at each sample (Figure 2, black).

Attenuation Correction
Beyond the discrete correction method published in (Richter et al., 2014), we developed an integral correction method to remove the above described attenuation effects.The procedure is based on the increase of the waveform intensity values in lower parts of the tree crown for portions of the pulse intensity, which have been reflected in higher parts of the crown.The attenuated received waveform can be corrected by stepwise increasing the amplitudes after each sample with an appropriate correction factor ci. Assuming that the attenuation of a laser pulse in the canopy fulfills the above described integral model, the correction factor can be developed by inverting the attenuation process.For this purpose, we calculate the proportion pi of the signal which was reflected at the first sample by setting the amplitude rwi of the received waveform at the sample point in relation to the reference value B ref (see figure 2, bottom): The attenuation correction factor ci is obtained as follows: In this way, the attenuated received waveform can be corrected by stepwise increasing the amplitudes with the appropriate correction factor ci.The deconvolution of emitted waveform and corrected received waveform results in the enhanced differential backscatter cross section, which is corrected for attenuation effects and can be inverted to the biophysical structure.Our correction method is based on the (in fact only confinedly valid) assumption that all surfaces, which interact with the emitted laser pulse, have a similar reflectance.Strictly speaking, this applies only on synthetic data because the reflectance of leaves, branches and forest floor differ from each other.The results in section 4, however, show that the quality of the volumetric forest reconstruction is nevertheless considerably improved.

Reference Selection
As described in section 3.3, we have to estimate the proportion pi of the signal which was reflected at the first sample to calculate the attenuation correction factor ci. To estimate the reflected proportion, information on the energy of the incoming laser pulse is required.Although the emitted waveform is recorded, a calculation of the incoming signal is not possible due to the unknown amplification of the receiver channels.Therefore, another reference value is needed.We assume that the maximal possible backscattered signal corresponds to the received waveform resulting from a propagation through the canopy without any interactions.The synthetic data set contains such a reference case (Figure 2) which make the reference value available.In case of real world data, a suitable value has to be derived from the recorded data set itself.For this purpose, we automatically detect all waveforms which fulfill the following requirements: • laser pulse penetrates to the ground without any canopy interactions (one peak on the ground) • waveform is located between vegetation waveforms (to exclude waveforms originating from forest tracks, clearings or fields) • waveform is recorded close to the nadir of the scan lines As can be seen in the histogram of all eligible reference waveforms in Figure 5 (top), the integrals considerably differ.One might either choose the mean value rvmean the areas under all eligible waveforms (for this data set 1516, σ = 251) or maximum value rvmax (here 2531) as a reference value.The reference value should actually be larger than the integrals of all other vegetation-waveforms because those are influenced by attenuation effects.The comparison with the histogram of all received waveforms in figure 5 (bottom) shows that this is not the case.The value derived in this way is more likely too small because the forest floor has mostly a lower reflectance than the leaves.Alternatively, the waveform with the highest integral (rvmax all = 2900), which results from the interaction with vegetation and forest floor, could be used.We have to take into consideration that the size of the correction factor significantly depends on the selection of the reference value.If it is too low, the correction factor becomes too large, and vice versa.In the process of validation we investigated the different variants for reference selection.The results are presented in section 4.

Validation with synthetic data
We generated the synthetic data as described in section 3.2.Figure 6 (top) shows the initial received waveform resulting from the convolution of synthetic emitted waveform and synthetic differential backscatter cross section as well as the received waveform after consideration the attenuation.The results of the attenuation correction are presented in figure 6 (bottom).Herein, initial waveform and the corrected waveform lie in top of each other, proving that they are identical (as expected).

Validation with real world data
In addition to the synthetic data, a real world data set recorded with a Riegl LMS-Q680 scanner in March 2010 in the area of Oberlausitz in Saxony (Germany) has been processed.We concentrate our investigations on a region covered with mixed forest.The comparison of the raw data and the corrected data after applying integral correction method shows the benefit of the developed correction method.As an example, the data of a section of one scan line is visualized as a vertical view with color coded waveforms (Figure 7).An improvement of the forest structure representation in the lower parts of the canopy can clearly be recognized.
Assessing the quality of the attenuation correction method, we have to take into consideration that the result of the correction significantly depends on the proper selection of the reference value B ref .Figure 8 shows the comparison of raw data and corrected data based on the different reference selection variants.For completeness, the correction result of the discrete correction method (Richter et al., 2014) with rvmax as reference value is also shown.
In general, the integral correction method provides a stronger correction than the discrete method.The data is only weak corrected with the discrete method.The chosen example highlights a major advantage of the integral method: the independence of peak detection.An important factor (and a weakness) in the discrete method is the reliability of the peak detection.The used minima/maxima search algorithm provides a simple and fast detection of peaks but is not able to find entirely every one.
Depending on the selected reference value, the correction results of the integral method vary in strength.Selecting rvmean, the correction turns out to be stronger than expected for many waveforms, even though this effect is not very distinctive in the chosen example.In some cases, the corrected amplitude of the ground reflection becomes higher than the amplitude of the reference waveform (that also result from ground reflections).Assuming a homogeneous reflectance of the forest floor, the amplitudes of both ground reflections should actually be equal.Obviously, the reference value is set too small in such a case.It is therefore recommended to prefer a larger reference value to avoid overcorrection.However, with rvmean all as reference value, the data remains rather weakly corrected for the majority of waveforms.As a compromise, the solution with rvmax seems most plausible.This statement has to be verified with further tests.

CONCLUSION
In this study, we developed an integral correction method to compensate for attenuation effects caused by reflections in higher regions of canopy.Generating volumetric forest stand representations waveform data leads to a significant improvement of the radiometric quality of voxel data.Thus, the accuracy of forest parameters derived from the volumetric representation will improve, especially for applications extracting information on the vertical structure of the lower canopy parts.Inevitably, the presented correction method is based on confinedly valid assumptions.Nevertheless, accepting those assumptions will usually lead to a better result than determining forest parameters from uncorrected data.
Further experimental investigations are needed to assess which correction approach conforms best to the physical reality.For this purpose, a recreated tree model of exactly known geometry will be scanned with a terrestrial full waveform laser scanner to evaluate the functionality of the correction methods.Additionally, better insight about the selection of the reference value can be derived from these measurements.

Figure 1 :
Figure 1: Cartesian voxel structure with projected differential backscatter cross section

Figure 5 :
Figure 5: Histogram of the integrals of all eligible reference waveforms (top), histogram of the integrals of all waveforms (bottom)

Figure 6 :
Figure 6: Initial and attenuated waveform (top), received waveform after applying the integral correction method (bottom)

Figure 8 :
Figure 8: Comparison of raw data and corrected data with different correction methods and reference values (rvmean/rvmax = mean/maximum value of the areas under all eligible waveforms, rvmax all = highest occurring integral considering all waveforms), green: attenuated waveform, red: corrected waveform, black: detected peaks

Figure 7 :
Figure 7: Vertical view of one scan line with color coded waveforms, top: raw data, bottom: corrected data (reference value: maximum integral of all eligible reference waveforms)