A VOXEL-BASED METHOD TO ESTIMATE NEAR-SURFACE AND ELEVATED FUEL FROM DENSE LIDAR POINT CLOUD FOR HAZARD REDUCTION BURNING

Drastic changes in the climate has revised the face of disaster management: it is contributing to abnormal intensity, frequency and duration of extreme weather and climate events. The year 2020 started with more than 100 fires burning across Australia. Hazard reduction burning has become a resolute and primary land management technique that contribute to the reduction of bushfire severity. One of the key variables to consider for this application is fuel load, as the accumulation of vegetation in a forest profile affects the intensity of the burn. Conventionally, fuel loads are measured by manually cutting the vegetation and physically measuring the quantity after dry heating. This process is expensive, and time consuming. There is an opportunity for these techniques to be digitised and automated to give results in a timely manner and work as a decision support tool for practitioners. This paper proposes a voxel-based approach that can be used for estimating fuel load and percentage cover of the vegetation, at the elevated and near-surface fuel/vegetation layer as a method to augment manual estimation. We use an airborne LiDAR pointcloud dataset of Vermont Place Park, Newcastle, Australia to test the method. The preliminary inspection of the results confirms the technique that can approximate conventional manual method. Next steps include performance testing including more dataset to derive quantitative measures on the approach.


INTRODUCTION
Australia has an extensive history with bushfire, with charcoal and fossil deposits indicating the presence of wildfire as far back as 400,000 years ago (Vigilante and Thornton 2016). Australia is known for its extreme weather events, however climate change is making these events increasingly severe (Steffen et al. 2019). The 2019/2020 bushfires have been dubbed as one of the worst bushfire seasons in Australia, burning more than 2 million hectares of land in New South Wales (NSW) and Queensland alone (Chang 2020). From 1967 until 2013, bushfires have cost the Australian economy 4.7 billion dollars (Handmer et al. 2018). Along with the loss of human lives and damage to the built environment, natural habitats of native species have been severely impacted, conditions made favourable for invasive species and the air quality is also affected, as illustrated in Figure  1. Climate condition, topographic features, fuel/vegetation condition and its accumulation are the factors that contribute to bushfire behaviour (Gould and Cruz 2012). Increasing temperature from global warming has prompted the worlds climate to change abruptly, making hot days hotter. Resulting in dryer vegetation, drier the vegetation the faster the fuel will ignite. Hazard reduction (HR) burning is a technique that is commonly used by fire practitioners worldwide, it is also known as controlled burning, prescribed burning, backfire, swailing or burn-off (Gage 2016). Its application is by deliberately reducing accumulated vegetation through direct burning in bushfire prone areas. In Australia, HR burns are conducted based on several guides, each state has its own guide to assess fuel. Commonly used guide in New South Wales (NSW) are the Overall Fuel Hazard Assessment Guide (OFHAG) and Knee-waist-shoulder Estimation Guideline (Hines, et al. 2010, Forestry Commission of Tasmania, 1984. OFHAG outlines a hazard rating from low to extreme for vegetation layers in a forest profile. One of the metrics used to assess hazard rating is Fuel Load (FL) -the accumulation of vegetation over an area (Volkova et al. 2016). FL is considered along with weather and topography when conducting hazard reduction burning (Penman et al. 2011). There are two conventional method for acquiring FL data, common attribute of both the techniques are that the fire practitioner will go to the dense forest profile and visually assess at what location has the highest accumulation of vegetation. The difference between the two is that: 1) Deliberately cuts the vegetation of interest, oven dries it and weigh the vegetation. 2) Fire practitioner goes to the site and assess the vegetation based on plant percentage cover and surface litter depth, at a certain radius around the assessor. Both of this technique are time consuming and costly. At the same time, this visually assessments are prone to errors due to the technique's subjectivity of the assessment (Volkova et al. 2016, Spits et al. 2017 Skowronski et al. 2011, Hermosilla et al. 2014, Chen, Zhu et al. 2016. Light Detection and Ranging (LiDAR) facilitate rapid acquisition of detailed data and allow tasks to be increasingly digitised and automated (Yebra et al. 2015, Price andGordon 2016). LiDAR has been widely used worldwide for classifying forestry characteristics. In Turkey (Inan et al. 2017) investigated the potential of airborne lidar to characterise the forest profile based on: tree height, tree diameter, crown diameter and the tree volume. This research quantified the forest stands (trees) in a forest profile. (Almeida et al. 2016) used a portable profiling LiDAR for fire susceptibility and contrasting fire damage in the central Amazon. To determine the individuality of trees from colour pointcloud data. However, both this research did not fully explore the potential of LiDAR. Gorte and Winterhalder (2004) and Thies et al (2004) reported a voxel-based approach using mathematical morphology, skeletonization, connected component labelling and shortest path computation to reconstruct 3D trees from terrestrial scans. Hermosilla, Ruiz et al. (2013) in USA, Spain and Canada developed a methodology for the estimating the forest structure and canopy fuel parameters. They developed metrices for input variables for predictive fire models using forward stepwise multiple linear regression. This study used small-footprint and large-footprint full-waveform LiDAR sensors that collected data for over a decade. They used two approaches to answer the limitation present from full-waveform data (unable to detect detailed structure of the forest) (Listopad et al. 2012): Pulse detection and waveform modelling (Mallet and Bretar 2009). Similar approaches were used in different ecological areas in USA to quantify forest structures and species (Riano et al.2004, Skowronski et al. 2011, Zhao et al. 2011, Gonza´lez-Olabarria et al. 2012. In Canada and France, Grau et al. (2017) estimated 3D vegetation data with terrestrial laser scanner data using voxels based approach. They used ray tracing method to simulate realistic TLS data acquisition, then simulated realistic vegetation scene from a 3D range of plant area index (PAI) distribution and calculated voxel area density. This research evaluated the potential of voxels, and lead to the development of a tool called VoxelSD. This tool can be used to assess 3D vegetation density from TLS scans. It is an innovative way to estimate the distribution of plant/ leaf area with accuracy.
Over 60 years of extensive research, Australia has demonstrated several systems that could model and predict fire behaviour (Gould and Cruz 2012): CSIRO Grassland fire spread meter (Cheney et al. 1998;Cruz et al. 2015b), Phoenix Rapid Fire (Tolhurst and Chong, 2008) and the oldest model being the McArthur (1967). One of the common variables needed for this system is the FL data. (Chen et al. 2016) developed two predictive models of forest surface fuel load based on: LiDAR indices, forest types and previous fire disturbances. They classified the forest fuel structural characteristics based on surface fuel depth and percentage cover at distinct layer using a Terrestrial Laser Scanner (TLS): they classified the fuel using a raster image of a surface fuel depth that is interpolated based on the TLS points height values (h) at surface fuel layer. They then generated a scatter diagram for forest vertical profiles by plotting the density of LiDAR points against height and used a bimodal distribution to identify LiDAR points. Two models were created: 1 st represents the density distribution of LiDAR points across vertical profile of understorey shrubs and the 2 nd component plots the density distribution in overstorey fuel. Locally weighted scatterplot smoothing (LOWESS) was applied to smoothen the scatter plot. Two models were compared and the cut point between the two components of the bimodal curve was utilized to stratify and characterise the vertical structure of the forest and derive LiDAR indices for different vegetation layers. This LiDAR derived stratification method provides a significant contribution in vegetation classification, forest habitat mapping as well as forest wildlife conservation. (Spits et al. 2017) investigated the surface and near-surface bushfire attributes from Image-based point clouds. The image-based point cloud was captured from smart phones loaded with Fuels3D application. Fuels3D application is a tool for measuring fuel hazard and fire severity in the forest understorey (op.cit).
Most studies done to characterise forest fuel are attributed to the canopy layer of the forest profile, understandable as its responsible for the rate of fire spread and severity (Gould, McCaw et al. 2011, Gould and Cruz 2012, Chen et al. 2016, Price and Gordon 2016, Spits et al. 2017. However, there are limited studies done in the understorey of the vegetation layer. In this paper, we propose a voxel-based approach to estimate the volume of fuel in near-surface and elevated fuel layers below the canopy. Voxelisation is a useful spatial analysis technique when applied to dynamic phenomena like wind, fire, air and noise pollution (Gorte and Zlatanova, 2016). This approach may assist in quantifying FL from pointcloud data by building a threedimensional image of the composition of the forest. The next section provides an overview of HR burning and explains the methods to quantify fuel load; then the paper elaborates on the voxel-based methodology to estimate the volume of surface and elevated FL. The paper concludes with proposed future research and development. Fuel can be classified based on five structural layers in the forest: canopy fuel, bark fuel, elevated fuel, surface fuel, near-surface fuel, as illustrated in Figure 2. At the canopy layer, leaves and twigs are usually found on the tallest layer of the forest or woodland. The canopy is included in the elevated fuel. Bark fuel is the flammable bark on tree trunks and upper branches. Elevated fuel shrubs and juvenile understorey plants are usually 2-3m in height. Near-surface fuel are grasses, low shrubs and heath, these often have suspended components like leaves, barks and twigs (Hines et al. 2010

Assessing Fuel Load
The conventional way to assess FL in NSW is based on the kneewaist-shoulder estimation from the HR burning under dry forests guideline (Forestry Commission of Tasmania, 1984) as illustrated in Figure 3. This guide estimates FL using a visual assessment of surface fuel. A fire practitioner will go to the site, locate a spot and visualise at 2m radius (or 10m radius) around the assessor and then estimate the surface litter depth of fuel that is <6mm in diameter. Once the surface litter is established it is then correlated with the percentage of plant cover in that 2m radius circle. The combination of the litter depth and plant coverage provides an estimate for FL.   The percentage coverage of plant matter is used to estimate FL. For every 10% of cover at 20mm surface litter depth approximates to 1 t/ha, for every 100% at 20mm surface litter approximates to 10 t/ha. If the plant coverage seems less than 100%, then the percentage cover needs to be lowered to 80% then the estimated FL would be 8 t/ha (Forestry Commission of Tasmania, 1984).

Fuel Hazard rating
The resulting FL estimates are converted to surface fuel hazard scores as low, moderate, high, very high and extreme for the area of interest. At coverage less than 10%, the hazard is considered low as it has no effect to fire behaviour. At 10% to 20% it is considered moderate hazard rating as it dictates the flame height. Greater than 20% onwards the hazard rating is very high and significantly contributes to the flame height and rate of spread of fire. Assessing fuel hazard rating at the elevated surface fuel level is slightly different from near-surface fuel. Below 20% fuel coverage it is considered low. From 20% to 30% it is considered moderate. At 30% to 50% is considered high. HR burning is required for elevated surface fuel after the plant coverage reaches 30%. The near-surface fuel is of high significance when conducting HR burning because when fire spreads, that layer will work as ladder for fire to travel to the other layers in the forest profile (Spits et al. 2017).
These concepts, dynamics and calculations are naturally volumetric and typically deal with small-scale measurements in relation to larger areas of bushland. The ability to gather data remotely and process these data volumetrically and quickly offers the capacity to augment and extend existing techniques and provides an avenue for rapid data collection and analysis to assist in mitigating potential disasters, and pre-planning HR burns that may be quickly executed given the constricted timeframes HR burning is possible.
To assess and manage bushfire risk, planning fuel treatment and controlling smoke plumes, it is crucial to measure the fuel load and its arrangement to inform a versatile range of fire management activities (Zhou et al. 1998, Bradstock et al. 2010, Gould and Cruz 2012, Duff et al. 2017).
Miscalculation of FL can lead to inaccurate data and can therefore have significant management implications.
Underestimating FL leads to erroneously eliminating areas in need of HR burning and potentially leaves stakeholders underprepared in the event of a wildfire. Overestimating leads burning where its unrequired; this has environmental, cost and logistical implications. Inconsistent and inaccurate estimation of fuel load could lead to unreliable data for smoke dispersion models, prediction of greenhouse gas emissions and fire line intensity. Accurate and consistent fuel hazard data is imperative for fire and land managers in their decision-making process for fire management practices.
This research focuses on automating the process of measuring FL for HR burning though a voxel-based approach using airborne LiDAR pointcloud data as input. The methodology was applied to a dataset representative of typical Australian bushland, obtained from Fire and Rescue NSW. The study area shown in Figure 5 is Vermont Place Park, Newcastle, NSW. The Airborne LiDAR pointcloud is collected by a low altitude drone platform and has an absolute accuracy of <50 mm RMSE at 50m range, with 3 returns, an RMS ranging error of 30mm and a scan rate of 420l points/s (1 return). Figure  6 illustrates the test data set.

Study area, dataset and software
The pointcloud was transformed into local coordinate system, by truncating the geographical coordinates with the minimum values for x,y,z. A voxel grid was then generated across the workspace. Processing in voxel space brings three major advantages to pointcloud space: first, the size of the data to be processed can be reduced significantly; second, the voxel space allows to establish connectivity horizontally and vertically, which facilitates analysis and third, the effect of different point density can be reduced, which will be favourable for fuel estimation. High point density from overlapping flights (Figure 7) may give the wrong indication of a higher vegetation density.    The voxel analyses are performed with in-house J-based software. J is a powerful general-purpose array-based language that is suitable to develop algorithms for exploring problems with matrices and higher dimensional arrays. J is suitable for this study as a voxel space is a typical example of large array. Alongside J, Cloud Compare is used to visually inspect the pointcloud data.

Voxel resolution
The resolution of the voxel space will influence whether voxels are empty or filled as seen in Figure 9. Since we are using the count of non-zero voxels to estimate the fuel load, the resolution of the voxel space is very important. If the resolution is too course, FL is overestimated, if the resolution is too fine, there is a processing and speed overhead. Figure 9: Significance of voxel resolution for the number of non-empty voxels. Figure 9 shows a 2D representation of a voxel space: A illustrates that all cells are non-empty; higher resolution as in B has 16 cells non-empty and one empty, C illustrates 21 empty cells if the resolution is four times finer. Given the FL assessment used a dimension of 0.50m to knee height, a voxel resolution of 0.40m was determined to be optimal in that it each structural fuel layer could be classified as belonging to a particular stratum of voxels and subdivided accordingly.

VOXEL SPACE ANALYSES
The workflow consists of the following five steps: 1. Terrain estimation 2. Canopy classification and removal 3. Tree trunk classification and removal 4. Near -surface and elevated fuel segregation 5. Estimation of percentage coverage for elevated and near-surface fuel to calculate hazard rating The dataset for this research has more than the required standard for enabling detailed mapping of the forest profile: however, the density of pointcloud is also variable across the forest profile. Therefore, a voxel-based approach provides an optimal solution to add consistency and operability across low to high point densities (Thies et al 2004).

Terrain
When generating the digital terrain model, a common problem associated with LiDAR scanning of large forest profiles is its inability to penetrate foliage, which occludes the laser beam. This results in incomplete coverages on the understorey of the forest profile. The quality of the digital terrain is critical and influences the following steps. The terrain is determined in the pointcloud through square grid, subdivided into tiles of 4m x 4m. The lowest point in those tiles are located and a Triangular Irregular Network (TIN) are created from it, at each tile. Larger tile has a higher probability to have terrain points and the triangulation will approximate the terrain shape. However, the triangulation may lead to omission of hills and gaps. Voxel space is then created above the normalised TIN.
The digital terrain model is then modified so that the ground level and normalised to a horizontal plane, enabling the structural layers to be easily separated according to their relative height above ground level. To normalise the height data, the lowest nonzero voxel in each vertical column is transformed to z=0 (Eusuf et al 2020). After the terrain is normalised, we proceed to classify the other components of the forest profile; canopy, tree trunks, elevated surface fuel and near-surface fuel.

Canopy classification
Canopy classification assumes the crowns of the trees form a high-density region of filled voxels in the voxelised point cloud.
The determination of this region is based on applying 3D mathematical morphology operations Dilation and Erosion in a mixture of object and background voxels, where an object is defined as collection of adjacent object voxels.
• Dilation makes objects larger by changing background pixels adjacent to the object into object pixels; holes inside objects and cracks between objects may therefore disappear, i.e. become part of the object. • Erosion makes objects smaller by turning object voxels near to the outer boundary of an object into background; holes become larger; cracks become wider, small or narrow objects may disappear.
Whether a background voxel is close to an object, resp. an object pixel is close to the boundary (i.e. the background) is determined by the Structuring Element (SE). In this study we use a simple SE that only considers nearest-neighbouring voxels. Starting from the voxelised point cloud (Figure 10) we first apply two iterations of dilations, which adds a layer of 0.80m around all objects in that space. As a result, holes of up to 1.60m diameter with the canopy region get completely filled up. Also, the top of the canopy becomes 0.80m higher, and its bottom 0.80m lower.
Other objects, such as trunks and understory objects will extend by 0.80m in all directions, but some space between these objects will remain. The terrain surfaces become 0.80m higher. Next, we apply four iterations of Erosion. This will make originally small objects, even though they had just been extended, completely disappear. However, the (now filled) block of canopy will shrink (by 1.60m from all sides), but the core remains. Finally, we apply again two iterations of Dilation. Now the canopy will grow back to its original size, but the objects that had disappeared before will not come back ( Figure 11).

Figure 12: Connected component labelling
To the result of the above we now apply Connected Component Labelling. This operation assigns a unique label (as the voxel value) to each collection of adjacent object voxels, i.e. to each object. After this, the size of each object (i.e. the number of voxels in it) can be easily established from a histogram, and only the very large objects remain as canopy ( Figure 12).

Tree trunk classification
The tree trunks naturally cut vertically through the voxel space, through the near-surface and elevated surface fuel and up to support the canopy. For this study, the trunks would have resulted in overestimation of data when calculating the voxel volume of near-surface and elevated fuel. An in-house routine identifies vertically contiguous voxel columns assuming they approximate tree trunks; if there is a large enough connected group of filled voxels, we expect this to be typical of a trunk. The routine looks for groups of 10 vertical voxels with at least 9 filled. If there is a gap, it is considered to be an artefact of occlusion. This set of voxels was then isolated from the dataset. Figure 13 shows the forest profile in a section of the voxels space having width of 10 voxels. The mid graphics shows the isolated tree trunks, and the lower shows just the remaining voxels. They can be considered vegetation that have to be assessed of fuel volume. The colour coding shows red to blue: red representing the more non-empty voxel in the section and blue as the less non-empty voxels.

Elevated surface and near-surface fuel
Having removed everything but the structural layers of elevated surface fuel and near-surface fuel, the two layers can be further analysed to estimate the amount of undergrowth of vegetation at each level. It is important to state that the lowest voxel in each column represents the ground plane-having no LiDAR points below the ground level. This is voxel zero. At a resolution of 0.40m, voxel 1 and voxel 2 represent near-surface fuel and elevated surface is upward from voxel 3. Figure 15: Elevated surface fuel at 0.40m resolution Figure 15 shows the elevated surface fuel visualised using CloudCompare at resolution 0.40m and Figure 16 shows the vegetation in the near-surface fuel level across the entire site. These data can now be analysed and plotted to reveal plan coverage density and distribution information.

Estimation of fuel coverage for elevated and nearsurface fuel
Figure 17: 0.40m resolution voxel plots, where voxels do/do not contain near-surface fuel (above) and elevated fuel (below). Figure 17 shows the extracted and isolated structural layers of concern to HR burn practitioners. The non-empty voxels containing are plotted in red and the white areas represent empty voxels (with no LiDAR points). The density and distribution of the fuel can be assessed even visually. Figure 18: Densities of fuel types in circular areas with 10m radius of near-surface fuel. Figure 18 and Figure 19 illustrate circular areas having a radius of 10m, similar to the NSW knee-hip-shoulder methodology and the corresponding coverage percentages: 0-20% (blue), 20-40% (green), 40-60% (orange) and more than 60% (red). Figure 19: Densities of fuel types in circular areas with 10m radius of elevated fuel.
The circular areas shown in the figures can readily be used by practitioners using OFHAG guide. For example. Blue areas indicate plan coverage less than 20%. Based on the amount of near-surface and elevated fuel, these areas do not require HR burning. The other areas represent plant coverage greater than 20%. These areas have to be checked and qualified for HR burning. At this point, a visiting fire practitioner can arrive at the site pre-prepared to ground truth the data and inspect the areas indicating the highest risk, and possibly direct any burns away from the blue zones.
To assess the research finding the test area was visited and the vegetation fuel was visually inspected. The observations confirmed the results obtained for elevated fuel. The images in Figure 20 illustrate vegetation in areas, which have high (red area in Figure 19) and low vegetation coverage (blue area nearby). The classification of near-surface fuel needs more investigation and it is too early to draw conclusions. Near-surface fuel is highly correlated to the elevated fuel and the tree canopy. If the canopy or elevated vegetation are dense, the near-surface vegetation is considerably less. Figure 20: Images from the test site for elevation fuel: more than 30% coverage (above) and less than 20% coverage (below).

CONCLUSIONS AND FUTURE WORK
UAV drones and LiDAR pointclouds have become practically ubiquitous, and a represent a bigdata resource. The use of voxels, as demonstrated in this paper has enabled a substantial value-add to increase the interpretability, (re)usability and practical application. The workflow developed is low-cost, processor efficient and fit for the purpose of serving as a decision-support tool to support fire practitioners by improving the currency of data and in developing an evidence-base for helping map and prioritise the large-scale operation of HR burning. It should be noted that this research is still work in progress, but the results are very promising. The next steps will focus on ground-truthing and calibrating the voxel analysis method, specifically improving the terrain estimates and near-surface fuel. Further advancements in the include incorporating further sources of information into the model; particularly topographic and slope information that had not been examined in this study.
Given the window of opportunity to execute HR burns is contracting dues to climatic changes, it is important fire practitioners have up-to-date and longitudinal intelligence to best plan and perform the complex operations at hand.