VOLUME ESTIMATION OF FUEL LOAD FOR HAZARD REDUCTION BURNING: FIRST RESULTS TO A VOXEL APPROACH

In 2019/20 over 100 severe bushfires burned across the continent of Australia. The severity of these fires was exacerbated by many factors, including macroclimatic effects of global warming and, at the meso and micro scales, land management practices. The bushfire phenomenon cannot be stopped, however better management practices can help counter the increasing severity of fires. Hazard reduction burning is a method where certain vegetation is deliberately burned under controlled circumstances to thin the fuel to reduce the severity of bushfires. Fuel load is an important parameter to assess when hazard reduction burning, as the accumulation of vegetation in a forest profile affects the intensity of the burn. Conventional methods of measuring fuel load are time consuming and costly, and therefore it becomes increasingly important to investigate automated approaches for assessing fuel loads. This paper provides an overview of hazard reduction burning while explaining the methods to quantify fuel load. Then the paper presents our voxel approach in estimating the volume of fuel loads. The first results regarding different voxel resolutions are reported and analysed. This paper concludes with future steps and developments. * Corresponding author


INTRODUCTION
Global warming has induced the world's climate to change relatively abruptly. This climate change has forced the world to transition into a new climatic state, where natural disasters are predicted to occur with increasing frequency and severity costing lives, destroying infrastructure and disrupting established cycle (Wahlquist 2019, Gomes 2020. For example, increasing temperature causes vegetation to dry at a faster rate, and the dryer the vegetation, the easier it is to ignite. This is resulting in bushfires with a high burn-severity (Kose, Nikos et al. 2008). Severe bushfires affected Australia in 2019/20, notably affecting the air quality in major cities (Figure 1). The fear of loss of life and property to bushfire is exclusively pronounced in wildland-urban interface (WUI), which are areas where human development integrates with bush and rural vegetation (Gould, McCaw et al. 2011). It is paramount that communities and fire authorities in bushfire prone regions of Australia understand fire related characteristics of nearby bush (forest, grassland, shrublands and other vegetation types, in order to comprehend the potential hazard of fuels and their related fire behaviour for more effective decision making. The study of bushfires includes landscape systems and time phases. It is a multi-stakeholder, multi-variable and multi-scale problem. Bushfire is a phenomenon that cannot be completely mitigated. However, its severity can be influenced. Over the last two decades, a range of methods have been developed to reduce the severity of bushfires. Hazard Reduction (HR) burning has become one of the resolute applications in the management of fire prone ecosystems worldwide, where certain vegetation is deliberately burned under controlled circumstances to thin the fuel and reduce the severity of bushfires (Vigilante and Thornton 2016). Common method of visual assessment to assess hazard score for HR burning is Overall Fuel Hazard Assessment Guide (OFHAG) (Hines, Tolhurst et al. 2010, Gould andCruz 2012). OFHAG assess vegetation as fuel and gives a hazard rating about the five structural layers in a forest profile (Hines, Tolhurst et al. 2010): canopy, bark fuel, elevated fuel, nearsurface fuel and surface fuel as illustrated in Figure 2.  (Hines, Tolhurst et al. 2010) These guides are used by fire practitioners, to develop maps of bushfire fuel hazard for fire risk analysis, prioritise HR operations, conduct post burn assessments to determine the effectiveness and provide inputs into fire behaviour prediction models (Volkova, Weiss Aparicio et al. 2019). However, according to Volkova, Sullivan et al. (2016) these were developed to assess fuel hazard and were not intended specifically for fuel load estimation, they just provide the correlation of fuel load to Hazard rating but are not accurate. Fuel Load (FL) is the accumulation of vegetation in a forest profile measured in tonnes per hectare (t/ha) (Gould and Cruz 2012). Accurate measurement of FL is needed when assessing a forest profile for bushfire hazard. The FL data is highly important in the application of HR burning as it correlates to hazard rating. For example, underestimating the FL could result in inaccurate estimates of low bushfire risk thus eliminating the areas that should be prioritised for HR burn, potentially leaving the residents and properties under prepared (Spits, Wallace et al. 2017). In contrast, overestimating the FL will lead to overestimation of fuel hazard that could have implications to inflating building cost in bushfire-prone areas (Volkova, Sullivan et al. 2016). Similarly, the estimates of greenhouse gas emissions from bushfire would be poorly predicted using such poor-quality data (Volkova, Sullivan et al. 2016). The conventional method of estimating FL for hazard rating is by physically removing vegetation and then sorting, drying and weighing the fuel. This method is time consuming, costly and requires high labour involvement. Therefore, it becomes increasingly important to investigate automatic approaches to acquire FL data. Rapid data acquisition via remote sensing technology such LiDAR allows many of these tasks to be automated (Yebra, Marselis et al. 2015, Price andGordon 2016). LiDAR has been extensively used in study of forestry characteristics. (Skowronski, Clark et al. 2011) scanned the vertical profile of a tree using LiDAR systems to estimate the FL of a canopy developed from field plots and allometric equations and use LiDAR datasets to predict the density and weight of the canopy layer of a tree. (Chen, Zhu et al. 2017) developed two predictive distribution models to accurately estimate forest surface fuel load with LiDAR data through multiple regression analysis. This study recognised an accurate and coherent method to estimate the spatial variations in forest surface fuel load. Grau, Durrieu et al. (2017) have utilised specific development methods to extrapolate information on vegetation density. Most studies related to quantifying vegetation are in relation to estimating the canopy layer. However, according to (Chen, Zhu et al. 2017, Grau, Durrieu et al. 2017 there are some limitations when capturing and scanning vegetation in a forest profile to capture and scan the structure of the vegetation in a forest profile: 1. Occlusion effect: When the beam is intercepted by canopy elements, the spaces behind those elements are not sampled. 2. Partial hit effect: the beam is intercepted by either a branch or a leaf, which causes multiple echoes for one emitted shot 3. Intensity uncertainty: the return of the beam is linked to several factors, optical properties, orientation of the component, distance from the scanner, noise, and instrument gain In this paper, we attempt to mitigate the above-mentioned problems associated with LiDAR by proposing a voxel-based approach. Voxel is a volumetric pixel, that provides spatial analysis for dynamic phenomena like wind, fire, air and noise pollution, visibility analysis or navigation (Gorte and Zlatanova 2016, Aleksandrov et al 2019. Voxels are expected to provide a solution to solving some of the issues linked with LiDAR data. Grau, Durrieu et al. (2017) developed a voxelisation method to assess the distribution of plant area from a Terrestrial Laser Scanner (TLS). They mapped a TLS point cloud in a voxel grid and computed the number of echoes inside the voxel. Quantifying near-surface and elevated fuel layers from point cloud data are complicated, as points in the structural layer could be missing due to the inherent limitations in LiDAR.
The paper is organised as follows: 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; finally, first results based on different voxel resolutions are reported and analysed; The paper concludes with future research and developments.

METHODS FOR BUSHFIRE HR BURNING
Fire has played a fundamental role in the evolution of Australia's biota and remains a key driver of many of its ecosystems. The early practices of Indigenous Australians for hunting and farming induced the Australian biota to be tolerant to fire; certain areas were deliberately burned for the purposes of hunting and farming (Bradstock, Hammill et al. 2010, Gammage 2011. Indigenous Australians practised a unique fire management technique known as firestick farming. Firestick farming is a practice to cool-burn certain vegetation deliberately to facilitate hunting and to change plant composition and animal habitat in a certain area. This kept the environment stable and fresh (Gammage 2011). HR burning is like firestick farming. Instead of hunting and cultivating vegetation, now it is used to reduce the severity of bushfires and for maintaining the balance in the ecosystem.

Factors of Bushfire
Fuel ignites when enough heat is given to a combustible source and combined with oxygen. Heat, fuel and oxygen combine to form a 'fire triangle' that keeps the fire burning (Nolan and Thornton 2016). The intensity and severity of a bushfire are depended upon temperature, wind speed, topography, fuel moisture and fuel load. From all these parameters the fuel load is the only aspect humans can influence. Therefore, this research aims to quantify the volume of fuel load from nearsurface and elevated fuel layers. To be able to quantify volume of fuel load, it is important to know the characteristics of fuel and classify based on these characteristics.

Fuel
The availability of forest fuel determines the amount of heat that can potentially be released in a bushfire. Fuel exists in numerous forms, sizes, states and arrangements. Fuels arranged with high level density with horizontal and vertical continuity promote the spread of flames (Hines, Tolhurst et al. 2010). They are usually: fine or coarse, dead or alive, woody or non-woody, surface or canopy (Gould and Cruz 2012). Table 1 outlines the characteristics of fuel from the five structural layers in a forest profile, the fuel included in the canopy layer are leaves and twigs and are usually found on the tallest layer of the forest or woodland. Bark fuel is bark on tree trunks and upper branches. Elevated fuel shrubs and juvenile understorey plants are usually 2-3m height. Canopy is normally less than 4m height and can be found included in the elevated fuel.

Fuel Classification:
Classifying fuel is critical for effective fire management as it provides a simple way to input extensive fuel characteristics into fire behaviour models and support various land and fire management needs (Gould and Cruz 2012). The main objective of fuel classification is to create and catalogue fuel attributes which ideally represent all possible fuel beds or fuel types for a region and their subsequent fire behaviour and effects. The amount of fuel present at the elevated and near-surface structural fuel layer suggests the behaviour of fire and predicts the rate of spread (Hines, Tolhurst et al. 2010, Gould, McCaw et al. 2011. Therefore, reducing the fuel at these layers could reduce the severity of bushfires. Elevated fuel is considered hazardous based on fuel continuity (horizontal and vertical), height, weight and ratio of dead materials and/or thickness of foliage. Very high elevated fuel hazard dictates the flame height and rate of fire spread. Fires may even spread in elevated fuels, even though surface fuel is wet. Near-surface fuels, however, are a lot different from elevated fuel because fire at this layer will spread no matter the quantity or weather conditions (Hines, Tolhurst et al. 2010). There is also limited study done on these layers as most research focuses on the canopy layer in the forest profile.
Analysing the components in this layer will help make better decisions when conducting HR burn. Hence, the motivation of this research is to quantify the volume of near-surface and elevated fuel load.

Fuel Load:
Fuel load is one of the factors that determine the severity of bushfires and hazard rating for HR burning (Hines, Tolhurst et al. 2010, Gould andCruz 2012). Calculating fuel load precisely can result in accurate estimates of fuel to assist in the planning of HR burn. The fuel load of surface fine fuel at low hazard is 4 t/ha, at moderate fuel load litter is 4-8 t/ha and at high is 12-20 t/ha. The cover of fuel determines how the fuel will travel. This process is not cost effective, time consuming and requires a lot of labour. This process becomes cumbersome especially when dealing with large-scale forest profile. Therefore, this research investigates an automated cost-effective and rapid approach to quantifying the volume of fuel load from a point cloud using a voxel grid.

Overall Fuel Hazard Assessment Guide
The Overall Fuel Hazard Assessment Guide (OFHAG) is a visual assessment guide for the Australian context, used mostly by Eastern states. It forms a basis for similar guides in other states. The assessment of fuels and their respective 'hazard' rating are commonly based on visual field assessments of different structural fuel layers. OFHAG aims to attribute hazard rating to vegetation based on the structure and continuity of fuels, live to dead ratio, height and size for each of the four fuel layers. Visual assessment of fuel provides a cost effective and rapid method to characterise fuels in individual forest layers as compared to direct measurement techniques (Watson, Penman et al. 2012, Spits, Wallace et al. 2017). However, visual assessments are subjective and can be vulnerable to inconsistency due to variability between assessors (Zhou, Robson et al. 1998, Sikkink and Keane 2008, Watson, Penman et al. 2012. For this research, the OFHAG will be used as a guideline to classify the fuels.

METHODOLOGY
This methodology describes the steps in developing an application for estimating the volume of FL from elevated and near-surface fuel layer. To compute the volume of near-surface and elevated fuel from point cloud data, we put the point cloud data through a voxel grid with specific resolution. Once the appropriate resolution is realised, the terrain of the study area needs to be classified to estimate the height from the ground to canopy to classify the height of each fuel in the structural layer of the forest. To calculate the height from the terrain to the canopy, the canopy layer needs to be classified first. After the height has been established, we classify the tree trunks of the study area. As we are interested in calculating the near-surface and elevated fuel load, the tree trunk and canopy needs to be excluded. The last step is calculating the volume of the filled voxels in order to estimate the volume of near-surface and elevated fuel load, assuming they provide a good estimate about the volume of near-surface and elevated fuel load.

Workflow
The first step of the workflow is to put the point cloud data in a voxel grid with appropiate resolution. The resolution is a crucial factor as it determines the size of each voxel in a voxel space.
The resolutions for this study are determined through visual observation, based on the least amount of 'noise'. The terrain is identified after thorougly visualising terrain representations at different resolutions. The resolution with the least amount of 'noise' will be considered appropiate for the study area. The height needs to be identified after classifying the canopy layer, to know at what particular height the near-surface and elevated fuel load are located in and elevated fuel load are located after classifying the canopy layer.

Figure 3: Detailed workflow
The forest profile needs to be horzontally sliced to be able to separate the five structural layers. Separating these layers will allow us to focus specifically to the layer that needs to be analysed; for this research its near-surface and elevated fuel layer. To prevent overestimation of point cloud data, the tree trunks in the forest profile need to be identified and classified. Once these steps are completed and the cut-point of understorey (surface fuel, near-surface fuel and elevated fuel) and overstorey (canopy & tree trunk) fuel layers are identified through horizontal slicing, we then estimate the elevated fuel and near-surface fuel volume by calculating the volume of filled voxels.

Voxels
A voxel is a 3D grid point in a closed axis-aligned unit cube in a Voronoi neighbourhood (Cohen-Or and Kaufman 1995)( Figure  4). A voxel is a volumetric pixel, a quantum of unit volume and has a numeric value integrated with it, which represents a property or independent variable of a real object or a value from a continuous field (Foley, Dam;et al. 1990et al. , Stoker 2009. A collection of voxels is described as a 3D array and represented as a voxel matrix (Stoker 2009).

Figure 4: Example of voxel in a voxel space
Voxelisation is a well suited approach for analysing point cloud data, due to these characteristics (Stoker 2009, Li et al 2018, Janecka, Karki et al. 2018) : • They are simple and unified • The grids can be thought of as a special graph with rectilinear connections  (Rosenfeld, 1981). Voxels can either be binary or numeric values (Stoker, 2009).
A 3D array of raster data records voxel values only, while the location of a voxel is defined implicitly by the position in the dataset (Gorte and Pfeifer 2004). This is in contrast to a vector data set, where the coordinates are recorded explicitly together with values and information, such as topology. In a 3D raster, 3dimensional phenomena are represented with different meaning attached to a voxel value. In the simplest raster representation of laser points from a point cloud, the voxel value range may be limited to 0; meaning the voxel is empty, and 1 meaning it contains one or more laser points. To test our methodology, a dataset describing typical Australian bushland was obtained from Fire and Rescue NSW. The study area is located at Vermont Place Park, Newcastle, Australia ( Figure 6). The obtained Airborne LiDAR point cloud has an absolute accuracy of < 50mm RMSE at 50m range, with 3 returns, RMS ranging error of 30 mm and a scan rate of 420k points/s (1 return).
The point clouds were analysed and voxelised through J programming (https://www.jsoftware.com/#/). J is a powerful general-purpose programming language, particularly suitable for developing algorithms for exploring problems where data are represented as matrices and higher-dimensional arrays. J is well suited for this study as it works efficiently with large arrays (millions to billions of rows). The LiDAR point cloud data we have has more than 50 million data points. The data points are represented as x,y,z with intensity. CloudCompare (www.cloudcompare.org) was the application used to visualise and analyse the point cloud data.

Assumptions
As mentioned above, the objective of this study is to identify fuel (combustible vegetation material) in the space between the terrain and the bottom of the tree crowns. The implicit assumption is that the presence or absence of material inside this elevation range is indicated by the presence/absence of points in the LiDAR point cloud. To this end we subdivide the entire 3D space of the study area in a 'voxel space', which is made up of a large collection of cubes of a uniform size. The space is a 3D rectangular block, of which the horizontal extent is defined by the study area (which is not rectangular, however) ( Figure 7). The extent of the voxel space in vertical direction ranges from just below the lowest terrain point in the area until just above the highest tree present.
Then, the LiDAR point cloud is transferred to the voxel space, after which cubes will be 'filled' or 'empty', depending on whether the LiDAR did or did not record a reflecting object in that part of the space. The above-mentioned assumption is that this corresponds to the presence or absence of material. Therefore, for example, many voxels in the tree crowns will be filled, because of the high density of leaves and branches. Also, voxels located at the terrain surface are expected to be filled, as real-world material is present there as well. Below the terrain, all voxels are expected to be empty, since no laser beams are reflected by underground material.
In the height region between terrain and crown bottom, however, we expect a mixture of filled and empty voxels, and the estimate of fuel load will be based on the ratio between their numbers.

Analysis of spatial resolution
An important parameter is the size of the cubes, i.e. the resolution of the voxel space ( Figure 8). We want voxels to be as small as possible in order to get accurate results, but when they get too small (so the number of voxels gets too large), many will remain empty simply because of the limited point density of the laser scanner making the measurement unreliable. The UAV-mounted LiDAR system that was used in the study records a point cloud with a density of (on average) 430 points/m 2 , which is very high in airborne LiDAR terms. Typical aerial survey point densities are around 20 pts/m 2 . An optimal voxel resolution can be established by looking at the voxels located at terrain surface: if there are (many) empty voxels at the surface, then apparently material is at a high risk of going undetected because the resolution chosen is too fine. We compared different voxel resolutions (Figure 9). At a resolution of 40cm we obtain a voxel space covering an area of 1,200 x 1,220 pixels of 40cm x 40cm = 234,240m 2 . However, only 127,000m 2 of this area is covered by the available data, consisting of almost 54 million points. We observe at this resolution that the terrain forms an almost closed surface in the voxel space ( Figure 9). This indicates that at almost every (x,y) position there are one or more laser beams reaching the 40cm x 40cm surface patch at that position, and there is very little chance that significant objects at or above the surface go undetected. Another consideration in this discussion should be how the fraction of laser beams penetrating the tree crowns (being available to detect combustible material in the understorey) depends on the chosen resolution, but this effect would have to be further investigated.

Classification of terrain points
Detecting the height of the terrain (i.e. constructing the DTM) is indeed the next step in the workflow: finding at every (x,y) in the voxel space the z of the lowest filled voxel. Small holes in the thus detected surface do occur as expected, and are filled up by morphological closing (Serra 1982). When analysing the voxel space column by column, i.e. one (x,y)-position at a time ( Figure 10 and Figure 11), and identifying the terrain as the lowest non-zero voxel in a column, we can normalize the voxel space by removing the (underground) voxels below that position, thus putting the terrain-surface voxel at position 0. After this, horizontal layers in the voxel space correspond to heights above the surface, irrespective of terrain relief ( Figure 12). After normalization the histogram of numbers of filled voxels per layer is meaningful, representing the vertical distribution of materials (re-sampled from 40cm to 1m bins in Figure 13).

Classification of canopy points
In a similar fashion as the terrain heights, the top of the canopy can be estimated: at every (x,y) position in the voxel space its height is given by the (z value of) the highest filled voxel. However, as opposed to the terrain, not every canopy voxel is expected to be filled. Therefore, we first perform a threedimensional morphological closing to the entire voxel space with a 5x5x5-voxel structuring element, in order to fill the gaps. This operation fills up holes of up to 160cm between filled voxels, thereby turning the entire crown-space into a connected block. After this both the top and the bottom of the tree crown region is delineated quite clearly ( Figure 14).

Classification of tree trunks
In the remaining understorey region, finally, we detect tree trunks, in order to separate them from the other material. Tree trunks are characterized by forming collections of (say, at least 10) filled voxels above each other (sharing the same (x,y) location). Because of possible occlusion we allow for maximum 20% of the voxels in such a collection to be missing (i.e. empty) ( Figure 15). Having identified filled voxels as either terrain, crown, or trunk voxels, the filled voxels that remain are designated to represent the combustible understory material we are looking for. The next step to be performed is calibration: linking this result to a fuel load estimate.

SUMMARY AND FUTURE WORK
This paper presents our voxel approach to estimate the volume of FL with appropriate resolution for the site. The tree trunks and the canopy were classified with a voxel resolution of 0.40m.
The classification of terrain and canopy was performed by investigating the voxel resolution and assumptions for lowest and highest possible LiDAR points in the data set. This process has been currently performed largely by human visual inspection. Further investigation & experiments are needed to confidently conclude on different voxel resolution depending on the density and type of vegetation as well as density of point cloud.
Classifying the tree trunks provided a promising result, as the tree trunk point cloud can overestimate the point clouds at nearsurface and elevated fuel layer. So, when estimating the volume of the filled voxel at elevated and surface fuel layer, and hopefully be able to estimate the FL at these layers. The next step in this research is classifying the remaining voxel space and computing the volumes of near surface and elevated vegetation.