Pre- and post-fire forest canopy height mapping in Southeast Australia through the integration of multi-temporal GEDI data, satellite images, and Convolution Neural Network

ABSTRACT

This study leveraged Convolutional Neural Network (CNN) models to estimate canopy height in Southeast Australian forests before and after the 2019–2020 bushfire event, using inputs from Sentinel-1, Sentinel-2, spectral indices, and terrain features and GEDI canopy height (rh99) as target variables. Our research consisted of three primary objectives: (1) an evaluation of GEDI rh99 canopy height in relation to an nDSM derived from airborne LiDAR, (2) the development of two CNN models for pre- and post-fire scenarios, and (3) utilizing the trained CNN models to generate pre- and post-fire canopy height maps and canopy height change maps across the continuous landscape. GEDI rh99 accuracy analysis revealed $R^2$ = 0.85, RMSE = 7.25 m and nRMSE = 0.23 for pre-fire GEDI and $R^2$ = 0.78, RMSE = 10.1 m and nRMSE = 0.3 for post-fire GEDI. The pre-fire CNN model achieved $R^2$ = 0.75, RMSE = 7.29 m and nRMSE = 0.22, and post-fire CNN model achieved $R^2$ = 0.62, RMSE = 11.67 m, and nRMSE = 0.34 when validated against nDSM. Bias analysis was conducted, revealing that GEDI rh99 underestimated nDSM across all canopy height ranges, and slope has no significant impact on GEDI accuracy. Notably, post-fire GEDI rh99 exhibited significant underestimation in high to extreme fire severity. The CNN models’ prediction displayed overestimation in short forests < 20 m and increasing underestimation in taller forests > 20 m and underestimation in high and extreme fire severity. Subsequent investigation indicated that canopy layer thinning in high fire severity conditions resulted in weak GEDI waveform signals and consequent bias. Additionally, the GEDI rh99 bias was propagated to the CNN models, together with insensitivity of spectral band values to canopy height in dense and burnt forests potentially contributed to models’ bias. Finally, we demonstrated CNN models’ ability to monitor forest recovery by generating canopy height maps spanning East Gippsland from 2019 to 2023.

KEYWORDS:

• Convolutional Neural Network (CNN)
• forest fire
• Global Ecosystem Dynamic Investigation (GEDI)
• Multi-temporal LiDAR
• Sentinel 1
• Sentinel 2
• forest height mapping

1. Introduction

Forest fires are a significant ecological challenge in Australia, leading to considerable impacts on human health (Yu et al. 2020), ecological disruption (Godfree et al. 2021), biomass and agricultural losses (Nolan et al. 2022), and increased carbon emissions (Ward et al. 2022). Assessing the impact of fires requires a comprehensive understanding of the forest structure both before and after the fire event, as these are closely linked to fire severity, forest recovery, and reburn fire severity (Gibson et al. 2020; McCarley et al. 2017; Skowronski, Gallagher, and Warner 2020; Wulder et al. 2009). Forests structure, namely forests with lower crown height, continuous vertical forest layers, and reduced canopy cover tend to experience the highest fire severity (Skowronski, Gallagher, and Warner 2020; Valentijn et al. 2019). The spatial distribution of pre-fire forest structure plays a critical role in identifying high fire-risk areas and predicting fire behaviour. Post-fire forest structure information is valuable to assess the fire impact (Laidlaw et al. 2022) and monitor the forest recovery (Meng et al. 2018; Viana-Soto et al. 2022) for effective fire management (Cansler et al. 2022).

Airborne LiDAR (Light Detection and Ranging) data has emerged as a valuable tool for acquiring high-resolution forest structure variables such as canopy height, canopy cover, tree density, and vertical biomass profile (Hermosilla et al. 2013; Kane et al. 2013; McCarley et al. 2017; Skowronski, Gallagher, and Warner 2020). When compared with post-fire, LiDAR data can accurately estimate fire-induced alterations in forest structure, such as canopy cover and vertical profile change, making it particularly useful for fire severity assessment (Hillman et al. 2021; Hu et al. 2019; Olga, Almeida, and Manuel Moreno 2020). While airborne LiDAR is a well-accepted tool for estimating forest structure alteration, access to multi-temporal airborne LiDAR data is limited, and where it is carried out, significant time gaps between successive surveys can occur. Limited access to airborne LiDAR data hinders the widespread applicability of accurate estimation of post-fire forest structure dynamics, which rapidly change and alter the structure-reflectance relationship over time (Jacobson and Jacobson 2010; van Leeuwen et al. 2010).

The spaceborne LiDAR data, Global Ecosystem Dynamics Investigation (GEDI), offers a solution to the poor spatial and temporal coverage of airborne LiDAR by providing repeated LiDAR measurements with global coverage. GEDI, stationed aboard the International Space Station since 2019, delivers extensive forest structure data through its full waveform LiDAR measurements (Dubayah et al. 2020). The GEDI Level 2A and Level 2B products offer structure metrics such as above-ground relative height, canopy cover, Plant Area Index (PAI), Plant Area Volume Density (PAVD) and Foliage Height Diversity (FHD), which enable canopy height estimation (Adam et al. 2020; Dorado-Roda et al. 2021; Liu, Cheng, and Chen 2021) and facilitate forest disturbance detection (Guerra-Hernández and Pascual 2021; Huettermann et al. 2023). When compared with airborne LiDAR, GEDI maintained good accuracy in canopy height estimation yielding an RMSE of canopy height less than 3.56 m in North American forests (Liu, Cheng, and Chen 2021) and an RMSE of 1.95 m − 3.96 m across the Mediterranean forests (Dorado-Roda et al. 2021). The GEDI data offers opportunities to observe the forest structure change and track forest recovery over time on a regional scale (Lin et al. 2023), which was normally not achievable by airborne LiDAR data.

Nevertheless, GEDI is scatter data and it has 25 m-diameter footprints, which means that it only samples a small area across the continuous landscape. Therefore, studies explored the feasibility of extrapolating GEDI-derived structure data using satellite imagery, e.g. Sentinel-1, 2, and Landsat, coupled with machine learning models such as support vector machine (SVM) (Gupta and Kant Sharma 2022; Morin et al. 2022), Random Forest (RF) (Myroniuk et al. 2023; Nandy, Srinet, and Padalia 2021; Potapov et al. 2021; Rishmawi, Huang, and Zhan 2021; L. Wang et al. 2020) and Deep Convolution Neural Network (CNN) (Lang et al. 2022). These studies achieved high mapping accuracy of forest structure over large areas. For example, global forest height estimations by combining GEDI data and Landsat imagery using an RF model were shown to achieve an RMSE of 9.07 m and $R^2$ of 0.61 (Potapov et al. 2021) and in another recent study global forest height estimation using GEDI and CNN (Lang et al. 2022) achieved an RMSE of 7.8 m.

Machine learning models capture complex relationships between spectral information and forest structure within GEDI footprints and project them over large areas. In particular, CNN models are characterized by the ability to extract hierarchical texture information from images, in addition to the spectral information, to solve complex non-linear relationships among the morphological characteristics of forests (Kattenborn et al. 2021). The CNN model captures the textural information and learns the relevant patterns relating to structure metrics. The ability of CNN models to map forest structure based on satellite image input has been shown to achieve good accuracy. For example, Lang et al. (2019) used CNN model to map landscape-continuous forest height in Switzerland and Gabon and achieved RMSE of 3.4 m and 5.6 m respectively. Becker et al. (2023) used CNN model to map Norway forest height and achieved RMSE and normalized RMSE of 2.298 m and 0.179 respectively. Thus, combining GEDI and other satellite inputs is a generally accepted method to extract forest structure along large spatial and temporal scales (Becker et al. 2023; Dirk Wegner 2019; Lang et al. 2022; and Lang, Schindler, Shah, Asif Manzoor, and Bais 2020).

While forest structure time frames have been generated to characterize global forest structure and height, few studies attempted to apply these approaches to map forest structure and corresponding structure change in the context of forest fires (Lin et al. 2023; Myroniuk et al. 2023) and provide validation result based on ground truth or more precise airborne LiDAR. This is perhaps because machine learning model requires extensive reference data for training, while bi-temporal (pre- and post- fire) LiDAR data is not commonly available, further compounded by the short-time window to measure the post-fire forest structure before it starts to recover. Therefore, this study aims to assess the feasibility of extrapolating pre- and post-fire GEDI-derived forest structure information, particularly canopy height, using a Convolutional Neural Network (CNN), in the context of the 2019–2020 ‘Black Summer’ fires in southeast Australia. Specifically, this study aims to

1. Assess GEDI-derived canopy height accuracy in Australian forests by comparing it to nDSM canopy height (nDSM) from airborne LiDAR.

2. Build pre- and post-fire CNN models to estimate pre- and post-fire canopy height and evaluate the model accuracy against nDSM canopy from airborne LiDAR.

3. Use the two CNN models to generate canopy height maps and canopy height change map in pre- and post-fire periods.

2. Study area

The study area encompassed the entire burn areas within Victoria and New South Wales states in Southeast Australia (depicted in Figure 1), spanning the geographic coordinates of 145.8° E to 154.1° E and −28° S to −38.1° S. This area was subject to devastating mega wildfires between November 2019 and February 2020, originating from the north and progressing southward. The burn area extent was determined by fire severity maps released by the departments of Victoria and New South Wales (Department of Environment, Land, Water & Planning, 2020; Department of Planning, Industry & Environment, 2020). These unprecedented fires impacted over 60,000 $\mathit{k}{m}^2$ of forest, with 25, 27, 28 and 20% areas experiencing low, medium, high and extreme fire severity respectively, leading to mean canopy height reductions of 2 m, 3 m, 4 m, and 9 m respectively. The study area is categorized as Major Vegetation Group (MVG) 3 (Department of Climate Change, Energy, the Environment, and Water, 2020), comprising predominantly of eucalyptus open forest, with canopy heights exceeding 30 m and canopy cover ranging from 30% to 70%, and an understorey of shrubs, ferns, and herbs. Additionally, study area in the north of New South Wales is characterized by a less humid climate, which is consisted of MVG 5, the eucalypt open woodland that features canopy heights below 15 m and canopy cover below 20%. These areas exhibited dense understorey shrubs, hummock or tussock grasses, and forbs.

Figure 1. Image (left) shows the burn area and fire severity classes based on fire severity maps of Victoria and New South Wales. Image (right) shows the airborne LiDAR survey areas.

3. Data preparation

3.1. Airborne LiDAR

Airborne LiDAR data was acquired from the Department of Environment, Land, Water, and Planning (DELWP) in Victoria. The survey spanned from April 2019 to November 2020, covering in East Gippsland of VICTORIA (Figure 1, right). The LiDAR scanner used for data collection was the Riegl VQ-780i (Riegel Laser Measurement Systems GmbH, Horn, Austria), which achieved a point density of 22.12 pts/$m^2$. The LiDAR data was divided into pre-fire (April – November 2019) and post-fire (February – May 2020) periods and normalized to above-ground height using the ‘lidar_ground_point_filter’ function from the ‘WhiteBoxTools’ Python package. Subsequently, the 99th percentile height was extracted from the normalized ALS data, and based on this value, a 10 m-resolution Normalized Digital Surface Model (nDSM) was generated to represent the reference canopy height.

3.2. Fire severity map

To investigate the GEDI and model accuracy in relation to known fire severity, we combined the fire severity maps of 2019–2020 from the natural resource management departments of Victoria and New South Wales (Department of Environment, Land, Water & Planning, 2020; Department of Planning, Industry & Environment, 2020 respectively). Both severity maps were based on similar modelling methods that utilized Random Forests models to classify post-fire Sentinel-2 images to the fire severity classes with reference data captured through the interpretation of high-resolution aerial images. The fire severity classes in the Victoria severity map were 2) unburnt, 3) low canopy scorch, 4) medium canopy scorch, 5) high canopy scorch and 6) canopy burnt, which corresponds to 2) unburnt, 3) low, 4) medium, 5) high and 6) extreme classes of New South Wales severity map.

3.3. GEDI data

GEDI L2A data covering the entire study area during the pre-fire (April – June, 2019) and post-fire (February – June, 2020) periods was collected. The GEDI L2A relative height metrics, rh99, which represents above-ground height at 99% percentile of the GEDI L1B’s waveform signal, was used to represent canopy height. GEDI data was filtered using the following rules: beam types = power, quality flag = 1, degrade = 0 and sensitivity > 0.95 following method in Dhargay et al. (2022). GEDI power beams were kept because they have full energy to penetrate deeper canopies (Dubayah et al. 2021) and subsequently have higher elevation and canopy height accuracy (Adam et al. 2020). Finally, a forest mask (Australian Bureau Of Agricultural And Resource Economics And Sciences 2018) and a burn area mask from the severity maps were used to filter GEDI points that fell within burnt forest areas.

3.4. Sentinel 1, Sentinel 2 and terrain maps

Sentinel-1 (S1) GRD, Sentinel-2 (S2) Level 2A, and terrain feature data were collected from Google Earth Engine (GEE). Collection dates were separated into pre-fire (March – June 2019) and post-fire (March – June 2020) periods.

The S1 GRD data from GEE underwent pre-processing using the Sentinel-1 Toolbox, including thermal noise removal, radiometric calibration, and terrain correction. Further processing involved border noise removal, speckle filtering, and terrain radiometric normalization (terrain flattening) based on the method from Mullissa et al. (2021), resulting in analysis-ready S1 data. The VV and VH bands were selected from the S1 data.

Regarding the S2 data, images with cloud cover over 30% were removed. Additionally, pixels of cloud probability surpassing 10% and cloud shadows were systematically removed from individual images. Spectral bands (B2, B3, B4, B8, B11, B12) denoting blue, green, red, near-infrared (NIR), and two short-wave infrared (SWIR) bands were chosen for analysis. Additionally, spectral indices, including the Normalized Difference Vegetation Index (NDVI), Normalized Burn Ratio (NBR), and Normalized Difference Water Index (NDWI), were computed from the S2 bands. The indices are known to be sensitive to vegetation dynamics (Pérez-Cabello, Montorio, and Borini Alves 2021; Veraverbeke et al. 2012) and are frequently used in height prediction models, fire severity assessment (Collins et al. 2018; Escuin, Navarro, and Fernández 2008; Kane et al. 2013), and post-fire forest structure modelling (McCarley et al. 2017).

In respect to terrain feature data, we used the 30-metre Digital Elevation Model (DEM) to generate a slope and aspect map. Our rationale behind incorporating slope data was to enable the model to discern the influence of slope-related bias on GEDI-derived canopy height (Adam et al. 2020; Liu, Cheng, and Chen 2021). Furthermore, the inclusion of aspect information aimed to empower the model to interpret spectral value disparities arising from the interplay between solar azimuth and slope aspect. Notably, terrain features have exhibited considerable predictability in modelling fire severity classes (Birch et al. 2015; Dillon et al. 2011), potentially enhancing estimation of post-fire canopy height in light of the relationship between heightened fire severity and foliage loss.

Finally, a median image was derived from each pre- and post-fire S1 and S2 image collection, and all image data was resampled to 10 m resolution with the Nearest Neighbour method, and a forest mask (Australian Bureau of Agricultural and Resource Economics and Sciences 2018) was applied to remove all non-forest areas. A comprehensive overview of input ALS, GEDI and image data is in Table 1.

Table 1. Summary of input image data, spectral indices, terrain features and GEDI data.

Image data	Bands
Sentinel-1	VV, VH − 10m
Sentinel-2	B2 (496.6 nm) − 10m
	B3 (560 nm) − 10m
	B4 (664.5 nm) − 10m
	B8 (835.1 nm) − 10m
	B11 (1613.7 nm) − 20m
	B12 (2202.4 nm) − 20m
Spectral indices	NDWI = (B8A–B11)/(B8A+B11)
	NBR = (B8A–B12)/(B8A+B12)
	NDVI = (B8A–B4)/(B8A+B4)
Terrain features	Slope, aspect − 30m
LiDAR data	Height metrics
GEDI	L2A rh99
Airborne LiDAR	99% percentile above-ground height

4. Method

4.1. GEDI canopy height accuracy assessment

The accuracy of GEDI generated heights was assessed by comparing GEDI’s rh99 with the airborne LiDAR-derived nDSM. To ensure the comparability of the two distinct datasets, median nDSM value within a 12.5 m radius (corresponding to the each 25-m GEDI footprint size) area were sampled at each GEDI sample point. Subsequently, the accuracy metrics, $R^2$ value, RMSE, and normalized RMSE (nRMSE) were computed between GEDI rh99 and nDSM. The nRMSE was determined by dividing the RMSE by the mean of nDSM, which yielded an adjusted RMSE across forests with varying heights.

A bias analysis was conducted by calculating the difference between GEDI rh99 and nDSM (hereafter referred to as the ‘GEDI rh99-nDSM difference’). This difference was grouped by the nDSM canopy height categories from 0 to 60 m with 5 m intervals and the slope categories from 0 to 30 ${ }^{\circ }$ with 5${ }^{\circ }$ intervals. For post-fire data, the GEDI-nDSM difference was grouped by fire severity classes, i.e. unburnt, low, medium, high, and extreme severity.

4.2. Training Convolution Neural Network

We developed a Convolutional Neural Network (CNN) model with five convolution layers to predict canopy height (as outlined in Table 2) based on satellite images. The model consists of four Conv2D layers, employing a 3 × 3 kernel size and 1 × 1 strides. The Rectified Linear Unit (ReLU) serves as the activation function for each convolution layer, effectively clipping negative values to determine relevant activated features. All convolution layers are padded to maintain consistent dimensions. Following every two convolution layers was a 2 × 2 max pooling layer to reduce feature map size and expedite training. The last pooling layer was followed by a 1 × 1 convolution layer, which increased the number of feature maps from 256 to 512. Global average pooling was then applied to transform the 2D output feature of the final convolution layer into 1D vectors, which was followed by fully connected layers for making the canopy height predictions.

Table 2. Summary of the CNN model. The output shape indicates the dimension of output feature maps.

Layer	Output Shape	# of parameters
Conv2D_1	(15,15,32)	3776
Batch_Normalization_1	(15,15,32)	128
Conv2D_2	(15,15,64)	18496
Batch_Normalization_2	(15,15,64)	256
Max_Pooling_1	(7,7,64)	0
Conv2D_3	(7,7,128)	73856
Batch_Normalization_3	(7,7,128)	512
Conv2D_4	(7,7,256)	295168
Batch_Normalization_4	(7,7,256)	1024
Max_Pooling_2	(3,3,256)	0
Conv2D_5	(3,3,512)	131584
Global_Average_Pooling	(512)	0
Dense	(1)	513
Total parameters	525,313
Trainable parameters	524,353
Non-trainable parameters	960

The CNN model was trained by minimizing the loss function including the L2 regularization to counter potential overfitting Equation (1).(1) $\mathit{Loss}=\frac{1}{n}{\sum }_{i=1}^n{\left(\widehat{y_i}-y_i\right)}^2+\mathit{\lambda }\times {\sum }_{i=1}^n{w}_i^2$(1)

(1) $\mathit{Loss}=\frac{1}{n}{\sum }_{i=1}^n{\left(\widehat{y_i}-y_i\right)}^2+\mathit{\lambda }\times {\sum }_{i=1}^n{w}_i^2$(1)

The part $\frac{1}{n}{\sum }_{i=1}^n{\left(\widehat{y_i}-y_i\right)}^2$ represents the loss function, mean squared error that takes predicted ($y_i$) and reference ($\widehat{y_i}$) heights as input. The part $\mathit{\lambda }\times {\sum }_{i=1}^n{w}_i^2$ is the L2 regularization, in which λ is the regularization parameter, which controls the strength of regularization and was set to 0.01, and $w_i$ represents the weights of the CNN’s layers. The L2 regularization helps in controlling the complexity of the model and can improve generalization to unseen data. Finally, for model compilation, the Adam optimizer was employed that incorporates adaptive learning rates and momentum to efficiently update the model’s parameters during training. An early stopping was called when the model has less than 1% accuracy improvement within 10 epochs.

To build the training dataset, S1, S2, terrain feature maps and three spectral indices were clipped by a 15$\times$15 pixels window (equivalent to 150 m × 150 m) at each GEDI point. The clipped images were then composed into tensors with a dimension equal to 15 × 15 x 13 (13 represents total image bands). This yielded a dataset of 900,000 tensors which were paired with GEDI rh99 and split into 70% (630,000) for training, 20% (180,000) for validation, and 10% (90,000) for testing. Both training and validation data were used in optimizing the CNN model, while the test dataset was utilized exclusively for evaluating the model’s performance. The preparation procedure was done for pre- and post-fire CNN models separately utilizing data from the corresponding periods.

The model accuracy was evaluated by the metrics $R^2$ value, RMSE and normalized RMSE (nRMSE) using the test dataset. Two accuracy results were generated by using GEDI rh99 and nDSM as validation height. In addition, bias analysis was conducted by calculating the difference between the predicted canopy height and nDSM (subsequently called CNN-nDSM difference). The difference was then grouped by nDSM height and slope categories same as the method in section 4.1 and fire severity classes. These factors are known well affected accuracy of image-based canopy height prediction model (Healey et al. 2020; Potapov et al. 2021; Zhanmang et al. 2020).

4.3. Canopy height mapping and bias distribution analysis

The trained pre- and post-fire CNN models were subsequently employed to extrapolate pre-and post-fire canopy height maps using satellite images. Images were divided into 15$\times$15 pixels subsets without intersecting the edge to avoid border effect and input into CNN models for canopy height prediction. The resulting canopy height maps and nDSM were then clipped to the Orbost_resurvey region (Figure 1, right, green), because both pre- and post-fire airborne LiDAR were available for this area, which allowed us to directly calculate canopy height change. Canopy height change maps from pre- and post-fire nDSM and CNN were generated to compare the distributions of canopy height gains and reductions. Furthermore, the difference between CNN height maps and nDSM was calculated to explore the distributions of prediction bias.

5. Result

5.1. GEDI canopy height accuracy assessment

Compared to the pre-fire nDSM, the pre-fire GEDI rh99 had accuracy of $R^2$ = 0.85, RMSE = 7.25 m and nRMSE = 0.23, and post-fire GEDI rh99 had accuracy of $R^2$ = 0.78, RMSE = 10.1 m and nRMSE = 0.3.

Bias analysis revealed that both pre- and post-fire GEDI rh99 overestimated nDSM by 2–4 m, in the canopy height range of 0–5 m (Figure 2(a)). For the pre-fire period, GEDI rh99 was mostly unbiased within 5–35 m range, but underestimation increased with height above 35 m. For post-fire period, GEDI rh99 showed a trend of underestimation within10–60 m, with a tendency of increasing underestimation towards taller forests.

Figure 2. Pre- and post-fire GEDI rh99-nDSM difference grouped by nDSM height range ($\mathrm{a}$), slope ($\mathrm{b}$) and fire severity ($\mathrm{c}$). GEDI-nDSM difference was calculated by subtracting GEDI rh99 by airborne LiDAR-derived nDSM.

There was no evidence to suggest that slope had a substantial impact on GEDI rh99 accuracy (Figure 2(b)). Pre-fire GEDI exhibited a slight underestimation of canopy height across all slope and height ranges, with the median of the GEDI rh99-nDSM difference approximately 1–2 m. Post-fire GEDI displayed an overall underestimation of 2–3 m, across all slope groups.

Unsurprisingly, there was no significant effect of fire severity on pre-fire GEDI accuracy (Figure 2(c)). Pre-fire GEDI rh99 exhibited 2–5 m underestimation in unburnt category and 0–3 m of underestimation in all other fire severity categories. In contrast, post-fire GEDI rh99 showed 2–3 m of overestimation in the unburnt category and an increasing underestimation trend toward higher fire severity. Strong underestimation was particularly evident in the high severity, with a median GEDI rh99-nDSM difference of 3.5 m, and the extreme severity category, with a median GEDI rh99-nDSM difference of 14 m.

5.2. CNN training result and bias analysis

The outcomes of pre- and post-fire CNN revealed distinct performance metrics. When validating against GEDI’s rh99, the pre-fire CNN achieved $R^2$ = 0.76, RMSE = 7.28 m, and nRMSE = 0.28. In contrast, the post-fire CNN model achieved lower accuracy, with an $R^2$ = 0.67, RMSE = 8.31 m, and nRMSE = 0.35. When validating against nDSM, the pre-fire CNN exhibited an $R^2$ = 0.75, RMSE = 7.29 m, and nRMSE = 0.22, while the post-fire CNN demonstrated an $R^2$ = 0.62, RMSE = 11.67 m, and nRMSE = 0.34. In general, the pre-fire CNN had significantly better accuracy compared to post-fire CNN, and post-fire CNN exhibited decreased performance when using nDSM as the validation height.

Bias analysis, revealed that both pre- and post-fire CNN tended to overestimate nDSM in short forests below 20 m and underestimate nDSM in tall forests above 20 m (Figure 3(a)). The most substantial underestimation was observed in the tallest height range of 55–60 m, where median CNN-nDSM difference reached 15 m for pre-fire CNN and 20 m for post-fire CNN.

Figure 3. Pre- and post-fire CNN-nDSM difference grouped by nDSM height range ($\mathrm{a}$), slope ($\mathrm{b}$) and fire severity ($\mathrm{c}$). CNN-nDSM difference was calculated by subtracting CNN predicted canopy height by airborne LiDAR-derived nDSM.

Slope did not have significant impact on the CNN model’s accuracy (Figure 3(b)). The pre-fire CNN showed underestimation around 3 m across all slope ranges, and the post-fire CNN demonstrated a constant underestimation of 8 m across 0–25${ }^{\circ }$ and 6 m in slope of 25–30${ }^{\circ }$. Fire severity had no effect on the pre-fire CNN (Figure 3(c)). However, post-fire CNN demonstrated an underestimation around 8 m in unburnt, low and medium fire severity classes, and this underestimation significantly increased to 10 m and 15 m in high, and extreme fire severity classes.

5.3. Pre- and post-fire canopy height, height change and bias mapping

CNN-predicted canopy height map demonstrated significant canopy height reduction across burnt areas from pre-fire to post-fire period (Figure ($\mathit{b},\mathit{d}$)), compared to minimal canopy height change observed in pre- and post-fire nDSM (Figure ($\mathit{a},\mathit{c}$)). Notably, logging areas identified by elongated shapes in pre- and post-fire nDSM (Figure ($\mathit{a},\mathit{c}$)) were not clearly depicted in the CNN-predicted canopy height maps (Figure ($\mathit{b},\mathit{d}$)).

Pre- and post-fire canopy height change map from nDSM demonstrated that significant canopy height reductions over 20 m were mainly concentrated in mountain valleys (Figure 4(e), shown as red dendritic patterns). In contrast, canopy height change map from CNN (Figure 4(f)) exhibited extensive canopy height reductions over 25 m, which coincided with areas of high-to-extreme fire severity on the right side of the map and were clearly distinct from areas of low-to-medium fire severity on the left side of the map.

Figure 4. Images showed pre- ($\mathrm{a}$), post-fire canopy height ($\mathrm{c}$) and height change ($\mathrm{e}$) from nDSM and pre- ($\mathrm{b}$), post-fire canopy height ($\mathrm{d}$) and height change from CNN prediction ($\mathrm{f}$).

Bias maps were produced by subtracting CNN-predicted canopy height maps by nDSM. In pre-fire period (Figure 5(a)), considerable canopy height overestimations over 15 m (indicated in blue) were prominently observed within the logging areas and in the middle of mountain valleys. These areas were found having canopy height approximately 10–20 m, aligning with the height intervals that had the most significant overestimation in the pre-fire CNN model (Figure 3). Underestimations (depicted in red) over 10 m were predominantly evident along the slope sides of the mountain valleys.

Figure 5. Images showed canopy height difference between CNN prediction and nDSM in pre-fire ($\mathrm{a}$) and post-fire ($\mathrm{b}$) periods.

Conversely, in the post-fire scenario (Figure 5(b)), distinct underestimations above 20 m were aligned with regions of high-to-extreme severity on the map’s right side, while overestimations below 8 m were predominantly observed on the map’s left side corresponding to low-to-medium fire severity and also seen in logging areas and mountain valleys.

6. Discussion

6.1. GEDI canopy height bias in relation to canopy height and slope

Our bias analysis indicated that both the pre- and post-fire GEDI rh99 underestimated airborne LiDAR-derived nDSM in short forests of 0–5 m, followed by a trend of increasing underestimation in forests over 5 m. Similar canopy height underestimation using GEDI rh100 was reported in the conterminous U.S.A. (C. Wang et al. 2022) for all height ranges. The study highlighted that the poorest GEDI accuracy was observed in broadleaf forests like Maple, Beech, and Birch, characterized by closed canopies and dense understory, that contributed to limited signal penetration and erroneous ground level detection. Such condition may explain the overall underestimation observed in the eucalyptus forests of our study area that were characterized by very dense canopies (median cover > 50%) and multi-layered shrubs in the understory.

However, in other studies of eucalypt forests, Dhargay et al. (2022) found a different pattern for GEDI accuracy results. In Victoria’s Central Highlands, Dhargay et al. (2002) reported that GEDI rh95 consistently overestimated canopy across all height ranges (0–60 m), and that the bias diminished from 5.79 m to 0.03 m in the tallest canopies. Similarly, Adam et al. (2020) observed GEDI rh100 overestimating canopy height in European temperate forests across all canopy height categories, with this overestimation decreasing from 5 m to 0.5 m for taller forests.

In addition, our GEDI bias analysis indicated that slope did not have significant impact on either pre- or post-fire GEDI accuracy, which is also not corroborated by findings of previous studies that indicated GEDI canopy height accuracy reduced in steeper slope (Adam et al. 2020; Dhargay et al. 2022).

6.2. GEDI canopy height bias regarding to fire severity

We also found that the post-fire GEDI rh99 significantly underestimated nDSM in the high and extreme fire severity classes but was unbiased in unburnt, low and medium fire severity. Upon close examination of the LiDAR point cloud and GEDI L1B waveform in extreme fire severity conditions, we found that the GEDI waveform signal was substantially diminished in forest layers that had been consumed by fire (Figure 6). It’s important to note that the GEDI L2A rh metrics are calculated based on the detected waveform, which signifies both the ground and canopy top levels using a signal thresholding method (Dubayah et al. 2021). In cases where forest layers were consumed and lacked sufficient reflected signal to surpass the threshold, the true canopy top level was not accurately detected, resulting in underestimation of canopy height. We identified that large underestimation was characterized by a unimodal GEDI waveform that reflected only the ground signals. This indicates that the number of detected signal peaks in GEDI waveform (corresponding to num_detectedmodes parameter in GEDI L2A) had a significant impact on GEDI elevation and canopy height accuracy, which aligned with the GEDI accuracy analysis in Wang et al. (2022). This phenomenon was frequently observed in more open forests where there were few tree trunks within a GEDI footprint, which happened very often after severe crown fire (high to extreme severity classes). Consequently, fire severity emerged as the primary determinant of canopy height underestimation in the post-fire GEDI data, causing lower accuracy than the pre-fire GEDI data. Changing the setting of the relative height algorithm (a1 to a6) in GEDI L2A can improve the accuracy, but the optimal algorithm varied across forest types (C. Wang et al. 2022).

Figure 6. Example of pre- and post-fire Airborne LiDAR point cloud ($\mathrm{a}$, $\mathrm{c}$) and GEDI L1B waveform data ($\mathrm{b}$,$\mathrm{d}$). The thinning of forest layer after fire led to diminished waveform signals, resulting in difficulty in detecting accurate canopy height.

6.3. CNN bias regarding to GEDI canopy height bias

The pre-fire CNN model demonstrated better overall accuracy than the post-fire CNN model. Notably, the pre-fire CNN displayed consistent results when using GEDI rh99 and nDSM as validation data, while the post-fire CNN accuracy dropped when validated against nDSM, indicating the challenges in modelling post-fire canopy height based on the GEDI canopy height. This challenge appears to stem from a bias within the original post-fire GEDI rh99 data that consistently underestimated nDSM values in high and extreme fire severity (discussed in Section 6.1). In subsequent experiments, we retrained the post-fire CNN excluding data of high and extreme fire severity. The result showed notable enhancement in performance, with $R^2$ = 0.67, RMSE = 9.82 m, and nRMSE = 0.27, comparing to the original post-fire CNN accuracy $R^2$ = 0.62, RMSE = 11.67 m, and nRMSE = 0.34. These findings provide compelling evidence of the impact from fire-severity-related bias of original post-fire GEDI to the model accuracy. Additionally, both pre- and post-fire CNN models consistently demonstrated a tendency to overestimate canopy height in shorter forests (<20 metres) and exhibited an increasing underestimation trend in taller forests, resembling the bias trend observed in GEDI rh99 (Figure 2). Consequently, the development of a method for accurate GEDI data filtering and GEDI canopy height calibration becomes important. Several prior studies have proposed GEDI data filtering methods (Lahssini et al. 2022; Lang et al. 2022; C. Wang et al. 2022). However, it’s essential to acknowledge that these filtering methods often reduced available GEDI data for model training. The influence of GEDI data filtering on model accuracy needs further investigation.

6.4. CNN bias regarding to limitation of image-based data

Our CNN models exhibited a prominent increasing underestimation of tall forests above 30 m. Such underestimation phenomenon was common in studies of canopy height modelling utilizing optical sensors, which was reported in the forests heights above 27 m in the United States of America (Rishmawi, Huang, and Zhan 2021), above 20 m in Sub-Saharan Africa (Hansen et al. 2016), above 40 m across continental Australia (Zhanmang et al. 2020), and above 30 m within global forest height estimations (Potapov et al. 2021). This tendency was often ascribed to the pixel values saturating above a canopy height level that linked to the relationship between image reflectance and the forest structure (Freitas, Mello, and Cruz 2005; Mutanga and Skidmore 2004). We found a constantly increasing or decreasing trend of S2 band value from 20 m to 40 m, which subsequently plateaued beyond 40 m, making spectral signature of forests above 40 m mostly similar. The pixel value saturation phenomenon is prominent in our study area characterizing tall eucalyptus forests (median height >30 m) with substantial canopy cover (median cover > 50%) and multi-layered shrub in understory layers, making image value insensitive to the alteration of canopy height.

Similarly, the backscatter values of S1 VV and VH bands exhibit saturation in forests exceeding 30 metres in height. This saturation is likely attributed to the limited penetration of S1’s short wavelength (C-band) through dense canopies. However, study utilizing PALSAR with longer wavelength (L-band), which offers higher penetrability, also found saturation within the height range of 20 to 30 m in Australian forests (Zhanmang et al. 2020). Our analysis of the correlation between PALSAR-2 HH and HV backscatters and GEDI-derived canopy height revealed that their relationships do not surpass those observed between Sentinel-1’s VV, VH backscatters and canopy height. We assume that the forest’s dense understory layers overall increase the volume scattering and makes it difficult to discern tall forests from the short forests. The issue of pixel value saturation, coupled with the inherent bias in GEDI rh99, has ultimately exacerbated the underestimation in CNN height predictions for tall forests.

Notably, strong underestimation shown in post-fire CNN in high and extreme fire severity was potentially attributed by the decorrelation between the spectral band value and canopy height. The pre- and post-fire ALS point clouds in a high-severity region showed that the top canopy height remained unchanged after fire (Figure 6) when the foliage of surface and elevated forest layers had been predominantly consumed. However, the three spectral indices, NBR, NDWI, and NDVI that feature the most distinct difference between burnt and unburnt forests, exhibited significant value drop from the pre-fire (median value equal to 0.74, 0.47 and 0.87) to post-fire (median value equal to −0.08, −0.12, 0.48) in the same area. The spectral indices value drop did not signify alterations in canopy height, instead reflecting the abundance of charred material and dead vegetation, leading to insensitivity of image value to actual canopy height. As a result, the tall burnt forests exhibited similarly low NBR, NDWI and NDVI to the low burnt forests and unburnt low forests that confounded the CNN model.

The model overestimation was widely observed in low forest areas such as logging areas and the centre of mountain valleys. Further investigation revealed that this bias was potentially caused by spectral signature similarities between low and tall forests. An example in a zoom-in area (Figure 7) shows the logging areas of elongated shapes displaying >15 m overestimation (Figure 7(a), green) which had canopy height around 10–20 m (Figure 7(b)) and had no distinguishable difference to nearby forests of height above 30 m in the NBR image (Figure 7(c)). These regions contained dense elevated forest layers that contributed to high values in indices like NDWI, NDVI, and NBR, generating a spectral signature similar to that of taller forests. This spectral resemblance was consistent across various bands in Sentinel-1 and Sentinel-2 imagery, compounding the challenge of effectively distinguishing between low and tall forests.

Figure 7. Plots shows ($\mathrm{a}$) pre-fire CNN-nDSM difference, ($\mathrm{b}$) nDSM and ($\mathrm{c}$) Normalized Burn Ratio (NBR). The elongated shapes were identified as logging areas, which were indistinguishable from nearby taller forests in NBR image.

6.5. Applications

Our CNN model integrating the multi-temporal GEDI canopy height and satellite images offers a valuable tool for generating time-series of canopy height across landscapes, which can greatly facilitate investigations into the dynamic interplay between changes in canopy structure and fire severity (Gibson et al. 2020; McCarley et al. 2017; Skowronski, Gallagher, and Warner 2020) and effectively trace the trajectory of post-fire forest recovery, which is crucial for research on forest recovery monitoring (Laidlaw et al. 2022; Meng et al. 2018; Viana-Soto et al. 2022). To illustrate this capability, we produced a series of canopy height maps for East Gippsland in Victoria, using satellite images captured during leaf-on season from 2019 to 2023 (Figure 8). These time-series maps clearly depicted prominent canopy height reduction from 2019 to 2020 in areas with high to extreme fire severity (Figure 8, bottom right) and followed by signs of recovery in 2021, 2022 and 2023. The CNN models can well facilitate studies requiring extensive spatial and temporal canopy height data, which can be challenging to obtain without access to airborne LiDAR data and field measurements, by solely relying on satellite imagery as input.

Figure 8. Canopy height maps predicted by CNN models using satellite images from 2019 to 2023 during March to June ($\mathrm{a}$ to $\mathrm{e}$) and fire severity map ($\mathrm{f}$) in East Gippsland. The year 2020 marked serious canopy height reduction (dark areas), which matched the high and extreme fire severity (bottom right), followed by a gradual recovery through 2021 to 2023.

6.6. Future works

The above discussion underscores the challenges arising from biases in GEDI data. and the inherent limitations of image-based data. Although the latter limitation is often unavoidable, enhancing the accuracy of GEDI data remains an area for improvement. Conventional methods of GEDI data selection based solely on quality, sensitivity, and degradation metrics may not be sufficient for deriving reliable canopy height data. Studies such as Lang et al. (2022) have introduced a comprehensive framework for estimating GEDI height uncertainty and identifying qualified GEDI data based on GEDI waveform (L1B) geometry. Additionally, research conducted by Wang et al. (2022) delved into GEDI canopy height accuracy in relation to factors like GEDI beam types, accuracy parameters, forest types, and time differences between GEDI and airborne LiDAR data collection, offering valuable guidelines for filtering qualified GEDI data using these criteria. Consequently, future works focus on selecting accurate GEDI data and calibrating GEDI-derived forest structure metrics by leveraging the characteristics of the GEDI waveform and quality-affecting factors to improve the CNN model accuracy.

7. Conclusion

In conclusion, this study comprehensively assessed GEDI canopy height accuracy in Southeast Australia, both before and after the 2019–2020 Black Summer fires. Pre-fire GEDI achieved an accuracy of $R^2$ = 0.85, RMSE = 7.25 m, and nRMSE = 0.23, while post-fire GEDI rh99 had an accuracy of $R^2$ = 0.78, RMSE = 10.1 m, and nRMSE = 0.3. Both GEDI datasets consistently underestimated canopy height across all categories, with a growing bias towards taller forests. Notably, post-fire GEDI data significantly underestimated heights in regions with high to extreme fire severity due to burnt canopy thinning, rendering it indistinguishable in the GEDI waveform.

We employed a Convolutional Neural Network (CNN) prediction model using satellite imagery (S1, S2, and terrain features) with GEDI rh99 as the target variable to estimate pre- and post-fire canopy height. The pre-fire CNN achieved an $R^2$ of 0.76, RMSE of 7.28 m, and nRMSE of 0.28 when validated against GEDI rh99, and an $R^2$ of 0.75, RMSE of 7.29 m, and nRMSE of 0.22 when validated against nDSM. The post-fire CNN achieved an $R^2$ of 0.67, RMSE of 8.31 m, and nRMSE of 0.35 when validated against GEDI rh99, and an $R^2$ of 0.62, RMSE of 11.67 m, and nRMSE of 0.34 when validated against nDSM. Both models overestimated heights below 20 m and exhibited increasing underestimation for heights exceeding 35 m. While slope had minimal impact on model accuracy, the post-fire CNN displayed significant underestimation in high and extreme fire severity regions. These model biases were potentially related to GEDI rh99 bias, as they displayed similar trends in canopy height, slope, and fire severity. Additional factors contributing to model bias included spectral saturation for forests over 40 m in height, spectral signature similarities between tall and short forests with dense elevated layers, and significant spectral changes unrelated to forest height in burnt forests.

Leveraging the trained CNN models, we generated time series forest height maps in regions with extreme fire severity and regional-scale height maps across Victoria and New South Wales states, effectively capturing the forest recovery process and spatially correlating areas of height reduction with regions of high fire severity. These findings highlight the potential utility of our models for forest recovery tracking and analysis of structure-severity relationships.

Acknowledgements

The first author gratefully acknowledges the Faculty of Science, Monash University, for providing the scholarships for the research. We would like specifically to acknowledge the Department of Environment, Land, Water & Planning (DELWP), the State of Victoria for providing the airborne LiDAR data in East Gippsland of Victoria and Alex Codoreanu from DELWP for the assistance in selecting and transferring the data.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The airborne LiDAR data used in this study are from Department of Environment, Land, Water & Planning (DELWP), State of Victoria in Australia. Permission for the other uses of the data needs to be sought from DELWP.