Forest Canopy Height Retrieval and Analysis Using Random Forest Model with Multi-Source Remote Sensing Integration

Abstract

Forest canopy height is an important indicator of the forest ecosystem, and an accurate assessment of forest canopy height on a large scale is of great significance for forest resource quantification and carbon sequestration. The retrieval of canopy height based on remote sensing provides a possibility for studying forest ecosystems. This study proposes a new method for estimating forest canopy height based on remote sensing. In this method, the GEDI satellite and ICESat-2 satellite, which are different types of space-borne lidar products, are used to cooperate with the Landsat 9 image and SRTM terrain data, respectively. Two forest canopy height-retrieval models based on multi-source remote sensing integration are obtained using a random forest regression (RFR) algorithm. The study, conducted at a forest site in the northeastern United States, synthesized various remote sensing data sets to produce a robust canopy height model. First, we extracted relative canopy height products, multispectral features, and topographic data from GEDI, ICESat-2, Landsat 9, and SRTM images, respectively. The importance of each variable was assessed, and the random forest algorithm was used to analyze each variable statistically. Then, the random forest regression algorithm was used to combine these variables and construct the forest canopy height model. Validation with airborne laser scanning (ALS) data shows that the GEDI and ICESat-2 models using a single data source achieve better accuracy than the Landsat 9 model. Notably, the combination of GEDI, Landsat 9, and SRTM data (R = 0.92, MAE = 1.91 m, RMSE = 2.78 m, and rRMSE = 12.64%) and a combination of ICESat-2, Landsat 9, and SRTM data (R = 0.89, MAE = 1.84 m, RMSE = 2.54 m, and rRMSE = 10.75%). Compared with the least accurate Landsat 9 model, R increased by 29.58%, 93.48%, MAE by 44.64%, 46.20%, RMSE by 42.80%, 49.40%, and the rRMSE was increased by 42.86% and 49.32%, respectively. These results fully evaluate and discuss the practical performance and benefits of multi-source data retrieval of forest canopy height by combining space-borne lidar data with Landsat 9 data, which is of great significance for understanding forest structure and dynamics. The study provides a reliable methodology for estimating forest canopy height and valuable insights into forest resource management and its contribution to global climate change.

1. Introduction

In recent years, the negative terrestrial ecological impacts of climate change have become more apparent and are very likely to intensify over the coming decades [1,2,3]. Forests constitute a crucial segment of land-based ecosystems and are instrumental in sustaining the cycle of carbon on a global scale. They play a key role in counteracting the effects of global warming, enhancing environmental conditions, and maintaining the equilibrium of natural ecosystems [4,5,6]. The measurement of canopy height in forests is a critical indicator of the overall structure of a forest and its capacity to store carbon. Large-area and high-resolution forest canopy height data are fundamental to regional and global forest carbon stock estimates. It is an important part of assessing the global carbon cycle and estimating forest carbon storage, biomass, and forest dynamics [7,8,9,10]. Traditional methods of measuring forest canopy height rely on field measurements [11]. Manual surveys based on resource inventories can provide accurate and high-resolution data, but they are resource-intensive, requiring extensive human, material, and time resources [10]. However, these methods are only applicable to small forest areas, posing challenges in obtaining forest canopy height data for large regions. Conversely, remote sensing technology has advanced, allowing for the rapid, accurate, continuous, dynamic, long-term, and non-destructive estimation of forest canopy height through the extraction of horizontal structural information from remote sensing data over large areas [12,13]. Various methods, including optical remote sensing, microwave remote sensing, and lidar, can be used for forest canopy height retrieval. These methods significantly improve the efficiency and accuracy of forest canopy height retrieval without causing damage to forest ecosystems [14].

The Landsat series of Earth observation satellites is the longest-running program in Earth observation, continuously monitoring the global land surface on a moderate scale. Millions of high-quality medium-resolution multispectral data have been acquired. These data are stored at receiving stations in the United States and around the world for research on global change [15]. The most recent addition to the Landsat program, Landsat 9, was launched in September 2021 by the National Aeronautics and Space Administration (NASA) and the United States Geological Survey (USGS) [16]. Landsat 9 is equipped with enhanced sensors, including the Thermal Infrared Sensor-2 (TIRS-2) and Operational Land Imager-2 (OLI-2). Unlike its predecessor, Landsat 8, the instrument provides radiometric and geometric data superior to those on the previous generation of Landsat satellites. Landsat 9 OLI-2 has a higher radiometric resolution, increasing the quantization from Landsat 8 from 12 bits to 14 bits, enabling the sensor to detect finer differences, especially in darker areas such as water or dense forests, and has higher imaging capabilities to capture horizontal structure data of large-scale forest canopy [17,18,19]. Compared to Landsat 8 thermal infrared sensors (TIRS), the TIRS-2 also significantly reduces stray light, resulting in improved atmospheric correction and more accurate surface temperature measurements. In addition, it contains abundant spectral features and vegetation indices, which can effectively reduce the resource consumption of field measurements. It is an important data source for estimating forest canopy height [16,19]. However, it faces challenges in acquiring vertical structural characteristics of forests due to poor penetration, saturation, low data quality, and substantial errors in estimating forest canopy height. Consequently, these factors impact the accuracy of the estimation [19,20]. However, space-borne lidar has a higher orbital altitude, wider observation range, and high penetration. These data provide a solution to the challenges faced by Landsat satellite image data. In addition, it also provides reliable data support for large-scale terrain elevation and forest canopy height inversion [21]. In 2018, NASA embarked on a groundbreaking endeavor by launching two laser altimeter missions: ICESat-2 (Ice, Cloud, and Land Elevation Satellite 2) and GEDI (Global Ecosystem Dynamics Investigation) [22]. GEDI, a cutting-edge multibeam waveform lidar, is tasked with measuring the global tree height and canopy density through its unique eight laser beam ground trajectories [23]. The system is the world’s first spaceborne lidar altimeter system for high-resolution measurement of three-dimensional vertical structure and topography of tropical and temperate forest vegetation [24,25,26]. ICESat-2 boasts the Advanced Topographic Laser Altimeter System (ATLAS), incorporating pioneering micropulse multibeam photon counting lidar technology, a first for satellite applications. The primary scientific objective is to gauge vegetation canopy height accurately [27]. ICESat-2 and GEDI have unique advantages and represent the most advanced level of the space-borne laser altimeter. They have significantly expanded our capabilities in understanding and studying forest biomass, carbon, and water cycles and have opened new avenues for biodiversity conservation [23,27,28].

Previous studies have successfully utilized ICESat-2 and GEDI data to estimate forest height [29,30,31]. However, due to the characteristics of multispectral data and space-borne lidar data, researchers have been exploring the combination and application of multiple sensors and measurements to retrieve forest parameters [32]. Space-borne lidar data directly provide precise information about the vertical vegetation structure [33]. When combined with optical remote sensing data or radar data, this combination enables the large-scale and comprehensive estimation of forest attributes, including height [20,34]. Qi et al. [35] explored the efficacy of simulating GEDI data for enhancing TDX InSAR data height estimation in temperate forests, mountainous coniferous forests, and tropical rainforests. The results demonstrated that combining both data types led to more accurate forest canopy height estimation than using just one data source [35]. Potapov et al. [36] and Tyukavina et al. [37] emphasized the significance of incorporating Landsat-based surface phenology indices into forest structure modeling. Potapov et al. [38] utilized a superpixel machine learning algorithm to generate a global forest canopy height map at a resolution of 30 m by integrating GEDI and Landsat 8 data. The findings indicated that, in comparison to ALS data, the RMSE and MAE were 9.07 m and 6.36 m, respectively. Wang et al. [39] integrated WorldView-2 stereo images with Landsat 7 data and employed six machine learning methods to estimate forest canopy height. The results demonstrated that the gradient-enhanced regression method provided the most accurate estimation of CHM, with the R of 0.64 and RMSE of 3.1 m. Additionally, the study highlighted the limitations of using solely very high-resolution stereo images for CHM estimation and suggested that integrating CHM from WorldView-2 stereo images with satellite vegetation indices could enhance the accuracy of CHM estimation. Ren et al. [40] demonstrated the accuracy of GEDI LiDAR data and multitemporal Sentinel-2 images in estimating forest diversity in large areas using four regression algorithms in machine learning. In a similar vein, Li et al. [9] integrated data from the ICESat-2 lidar with observations from Sentinel-1, Sentinel-2, and Landsat-8 satellites. They employed deep learning techniques and random forest models to generate a high-resolution CHM. Furthermore, Xi et al. [41] utilized a combination of ICESat-2 lidar data and Sentinel-1 and Sentinel-2 data. They employed random forest and gradient-boosting decision tree methods in machine learning to develop a forest canopy height estimation model. Likewise, Guan et al. [42] employed GEDI and Landsat 8 data to generate canopy height maps for estimating forest age. Zhu et al. [43], based on the combination of GED lidar data and Landsat 8 and Landsat 9 data, developed a forest canopy height estimation model using the BP neural network in machine learning. The results show that the Landsat 9 model is superior to the Landsat 8 model, and the combination model of GEDI and Landsat 9 is the best. These studies collectively demonstrate the feasibility and effectiveness of combining space-borne lidar data with multispectral data for forest parameter retrieval, using machine learning methods to achieve significant results. Despite the advancements in space-borne lidar technology represented by GEDI and ICESat-2, research on the synergistic application of the most recent Landsat 9 data with these lidar systems is limited. Although the combination of GEDI and Landsat 9 data has been discussed in the literature [43], the time consistency processing of Landsat 9 images in the literature still needs to be optimized. At the same time, the extraction of Landsat 9 data feature variables in the literature only analyzes part of vegetation index variables, and there is no more comprehensive selection and analysis of Landsat 9 data feature variables. To date, there have been few studies on the potential integration of ICESat-2 and Landsat 9 data to assess forest canopy height. GEDI and ICESat-2 stand at the forefront of space-borne lidar, yet the real-world performance and benefits of combining these data with Landsat 9 for detailed forest canopy height retrieval have not been fully evaluated or discussed. These gaps highlight the great opportunity to enhance forest canopy height research and improve understanding of forest structure and dynamics by optimizing data processing, assessment, and analysis, integrating advanced data sources, and using superior machine learning methods.

Therefore, this study employs GEDI data, ICESat-2 data, Landsat 9 OLI-2 data, and SRTM terrain factors as the primary data sources. Various combinations of data sources are employed and compared to individual data sources. The random forest regression algorithm is employed to model the retrieval of forest canopy height within the study area. The canopy height model (CHM) extracted from airborne lidar is integrated as reference data. Lastly, the accuracy and precision of models using single and multiple data sources are compared and analyzed. In this study, the superiority and reliability of forest canopy height retrieval results of the random forest regression model based on multi-source data integration were fully evaluated by combining Landsat 9 data with different models of space-borne lidar data. It is proved that it is possible to obtain forest canopy height with low cost, high accuracy, and high efficiency in a large area, which provides a reference for regional scale forest remote sensing monitoring.

2. Materials and Methods

2.1. Study Area

This research concentrates on the two forested areas located in the northeastern United States. The chosen study areas were Harvard Forest (HARV) in Worcester County, Massachusetts (42°32′24″ N, 72°10′12″ W) (see Figure 1). HARV serves as a representative rural/wild region in the Northeast, with elevations ranging from 160 to 415 m and an average elevation of 348 m. Its climate is classified as temperate–continental–humid. The main vegetation in this region is an Eastern deciduous temperate forest, primarily featuring red maple (Acer rubrum) and white pine (Pinus strobus) [44,45]. This study involved data collection from the selected study areas, as well as the development and testing of relevant methodologies.

2.2. Research Data

2.2.1. NEON Airborne LiDAR Data

The airborne lidar data used in this study were obtained from the National Ecological Observatory Network (NEON). NEON’s airborne observation platform is equipped with an imaging spectrometer, an airborne lidar, and a high-resolution color camera. This initiative aims to gather a range of data, including multispectral images, lidar measurements, and varied ecological indicators, to investigate the impacts of climatic shifts, alterations in land utilization, and the proliferation of invasive species on ecosystems at a continental scale [46,47,48,49,50]. The airborne lidar-specific parameters can be found in

Table 1. Specifications of the NEON airborne lidar instrument and its derived products [22].

NEON Airborne Lidar Instrument and Product Specifications
Laser wavelength	1064 nm (near IR)	Elevation accuracy	5–35 cm
Laser power	250 µJ	Derived products	DTM and CHM (in the format of 1 km by 1 km mosaic tiles)
Laser repetition rate	33–167 kHz	Product resolution	uniform grid (1 m × 1 m)
footprint diameter	0.25 m (at 1000 m flying height), 0.8 m in wide beam-divergence mode	Terrain parameters	elevation and slope
Sampling density	1–4 points per square meter	Vertical datum	GEOID12A
Horizontal accuracy	5–15 cm	Canopy parameters	canopy top height, relative height (RH), and canopy cover

[22]. The NEON lidar-gridded data products used in this study were primarily acquired in August 2022 and can be downloaded from the following address: Explore Field Sites|NSF NEON|Open Data to Understand our Ecosystems (neonscience.org) [51]. The canopy heights for reference forests were determined by selecting the 90th percentile value within the 25 m × 25 m grid of the canopy height model (CHM) [36].

2.2.2. GEDI Data

Table 2. GEDI parameters name and size [22,24,52].

Parameters Name	Parameters Size	Parameters Name	Parameters Size
Track height	~400 km	Footprint	25 m
Coverage	51.6° N~51.6° S	Geolocation error	8 m
repetition rate	242 Hz	Along-track distances	60 m
Pulse width	15 ns	Across-track distances	600 m
Wavelength	1064 nm	Product tested in this study	L2A (elevation and canopy heights)

[22,24,52] provides a record of the parameter names and sizes used in the GEDI system. The GEDI system consists of three lasers. Two of the lasers operate independently at full power, while the third laser is bifurcated into two beams, generating four beams. These beams are then used to generate eight light spot trajectories on the surface of the Earth through beam jittering. The radar records the complete waveforms of the footprints [24], which can be seen in Figure 2. For temporal consistency with airborne LiDAR data, this study downloaded the latest GEDI L2A version 2 products from January 2022 to December 2022 via the NASA Earth Data website (https://search.earthdata.nasa.gov/) (accessed on 12 October 2023) [53,54].

We filter out invalid GEDI footprint data by utilizing parameters inherent in GEDI L2A data to ensure high-quality footprint points are selected for training and validation. From the GEDI L2A product, we filtered the GEDI footprint within the study area based on latitude and longitude. This filtering process is guided by the following parameters: quality_flag = 1, rx_assess_flag = 0, rx_algrunflag_flag = 1, degrade_flag = 0, and sensitivity ≥0.9 as a method of selecting the desired footprint [52]. The optimal algorithm was analyzed to obtain the final valid GEDI dataset by evaluating the data for the six algorithms (see

Table 3. The parameter configurations for GEDI L2A’s six algorithms with σ being utilized to denote the standard deviation of the background noise level [29,31].

Algorithm Setting Group	Smoothing Width (Noise)	Smoothing Width (Signal)	Waveform Signal Start Threshold	Waveform Signal End Threshold
a1	6.5σ	6.5σ	3σ	6σ
a2	6.5σ	3.5σ	3σ	3σ
a3	6.5σ	3.5σ	3σ	6σ
a4	6.5σ	6.5σ	6σ	6σ
a5	6.5σ	3.5σ	3σ	2σ
a6	6.5σ	3.5σ	3σ	4σ

) within the GEDI L2A footprint and determining the common valid footprint among these algorithms [29,31,54].

2.2.3. ICESat-2 Data

ICESat-2/ATLAS includes 21 standard data products ATL01~ATL21, which are divided into five levels, represented by Level 0, Level 1, Level 2, Level 3A, and Level 3B, respectively [27,55]. Among them, the data products related to forest vegetation include ATL03 and ATL08. ATL03 provides a global positioning photonic data product, and ATL08 is a terrestrial vegetation data product that provides topographic elevation, canopy height, and relative height indicators every 100 m along the route through the official processing of ATL03 data. In addition, ATL08 can be associated with ATL03 products to provide the opportunity to calculate ground elevation, canopy height, and relative height indicators at different scales [41,56,57,58,59].

ATL03 and ATL08 version 6 data products are used in this study. Considering the acquisition time of airborne lidar data, GEDI, and Landsat 9, the inaccuracy of forest canopy height inversion caused by time difference is avoided. The data acquisition time is from 1 January 2022 to 31 December 2022. The data come from the National Snow and Ice Data Center of the United States (https://nsidc.org/data/atl08/versions/6) (accessed on 12 October 2023). Details are shown in

Table 4. ICESat-2 data acquisition time.

Data Type	Acquisition Time	Number of Documents
ICESat-2 ATL08	2022/01/01–2022/12/31	4
ICESat-2 ATL03	2022/01/01–2022/12/31	4

. To facilitate subsequent analysis, we extracted a set of relative height indicators ATL08 data: RH25, RH50, RH60, RH75, RH85, RH90, RH98, and RH100. For ATL03 data, we obtained signal photons by matching them to corresponding ATL08 data products and extracted forest height metrics from these signal photons for each 25 m swath.

To ensure the consistent analysis of elevation data from various datasets, we utilized NOAA’s VDatum, which can be accessed at https://vdatum.noaa.gov/vdatumweb/ (accessed on 12 October 2023), to align the GEDI data and ICESat-2 data vertical datum with that of NEON. It is worth noting that crown height remains relative to the ground and is unaffected by the chosen vertical datum [22].

2.2.4. Landsat 9 Data

Because discrete GEDI and ICESat-2 footprints do not provide continuous coverage, it is necessary to prepare optical remote sensing data for omnidirectional forest height estimation. Mid-resolution Landsat 9 data are improved based on the previous generation. They have the advantages of availability, wide coverage, and high accuracy, but they have not been widely used in the study of omni-directional forest height estimation. This study mainly uses Landsat 9 optical remote sensing data. A total of 73 Landsat9 Level 1T data images covering the study area from January to December 2022 were acquired from the Google Earth Engine (GEE) website platform (http://developers.google.cn/) (accessed on 12 October 2023).

The GEE cloud removal algorithm was applied to the images, and an annual synthesis of time series images was performed to generate Landsat 9 data covering the study area in 2022 through median processing [60] (see Figure 3). Based on spectral characteristics, texture characteristics, and the vegetation index, the canopy spectral reflectance information is provided using seven bands of coastal, blue, green, red, NIR, SWIR1, and SWIR2 after image processing of Landsat 9 data (see

Table 5. Landsat 9 operational land imager-2 parameters [61].

Band No.	Band	Band Range/µm	Spatial Resolution/m
B1	Coastal	0.43~0.45	30
B2	Blue	0.45~0.51	30
B3	Gree	0.53~0.59	30
B4	Red	0.64~0.67	30
B5	NIR	0.85~0.88	30
B6	SWIR1	1.57~1.65	30
B7	SWIR2	2.11~2.29	30
B8	Pan	0.50~0.68	15
B9	Cirrus	1.36~1.39	30

) [61]. Six kinds of feature variables are extracted: the single-band factor, multi-band factor, vegetation index, principal component analysis (PCA) factor, minimum noise fraction (MNF) transform, and texture factor [62,63,64,65]. Further details of all Landsat 9-related variables extracted in this study are summarized in

Table 6. Summary of characteristic variables calculated from GEDI L2A, ICESat-2, Landsat 9, and SRTM data [62,63,64,65].

Type	Characteristic Variable	Description
GEDI L2A	RH25, RH50, RH60, RH75, RH85, RH90, RH98, and RH100	GEDI extracted relative elevation (25th, 50th, 60th,75th, 85th, 90th, 98th, and 100th)
ICESat-2	RH25, RH50, RH60, RH75, RH85, RH90, RH98, and RH100	ICESat-2 extracted relative elevation (25th, 50th, 60th, 75th, 85th, 90th, 98th, and 100th)
Landsat 9	B1, B2, B3, B4, B5, B6, and B7	Landsat 9 bands 1, 2, 3, 4, 5, 6, and 7
	B24	$B2/B4$
	B74	$B7/B4$
	B76	$B7/B6$
	B345	$B3\cdot B4/B5$
	EVI	$2.5\left(\frac{B5-B4}{B5+6B4-7.5B2+1}\right)$
	DVI	$B5-B4$
	SLAVI	$\frac{B5}{B4+B7}$
	VI3	$\left(B5-B6\right)/\left(B5+B6\right)$
	PVI	$\sqrt{{\left(0.355B5-0.149B4\right)}^2+{\left(0.355B4-0.852B5\right)}^2}$
	NDVI	$\left(B5-B4\right)/\left(B5+B4\right)$
	RDVI	$\raisebox{1ex}{$B5-B4$}\!\left/ \!\raisebox{-1ex}{$\sqrt{B5+B4}$}\right.$
	ND43	$\left(B4-B3\right)\cdot \left(B4+B3\right)$
	ND67	$\left(B6-B7\right)\cdot \left(B6+B7\right)$
	PC1, PC2, and PC3	The first to third bands of principal component analysis.
	TCB	The tassel cap transforms the Brightness band.
	TCG	The tassel cap transforms the green band.
	TCW	The tassel cap transforms the wetness band.
	MNF1, MNF2, MNF3, and MNF4	Minimum noise fraction first to fourth band.
SRTM	elevation	Elevation extracted from DEM.
	slope	Slope extracted from DEM.
	aspect	Aspect extracted from DEM.

.

2.2.5. Auxiliary Data

Due to factors such as sunlight, soil moisture, and topographical changes, plant growth varies at different slopes and orientations, resulting in different canopy heights. To mitigate the effects of slope, aspect, and elevation on terrain and canopy height, we introduced ancillary data called the Shuttle Radar Terrain Mission (SRTM) [66,67]. SRTM Digital Elevation Model (DEM) data are available in two resolutions: 30 m and 90 m [68]. Topographic parameters such as elevation, slope, and aspect were extracted from SRTM data with a resolution of 30 m in 2022 obtained from http://earthexplorer.usgs.gov/ (accessed on 12 October 2023). To ensure the consistency of extracted variables, this study uses the central coordinates of GEDI and ICESat-2 footprints as the reference for variable extraction, effectively solving the problem of resolution differences among SRTM data, GEDI data, ICESat-2 data, and neon data sources [69]. Taking these parameters as input variables of the random forest model, the effect of terrain on forest height is effectively mitigated.

2.3. Methods

In this study, the GEDI, ICESat-2, Landsat 9, and SRTM data were used to extract the characteristic index (see Table 6). Then, the multivariate importance analysis is carried out on these characteristic indexes to obtain the final characteristic variable set (see

Table 7. The feature variables utilized in the final experimental modeling process [62,63,64,65].

Model	Characteristic Variable	Number of Characteristic Variables
GEDI	rh98, rh100, rh90, rh85, rh25, and rh75	6
Landsat 9	MNF2, B3, TCW, B74, B7, MNF3, B6, ND43, MNF4, B76, EVI, B1, B24, ND67, B2, B345, B4, pc3, and VI3	19
GEDI and Landsat 9	rh98, rh100, MNF2, rh85, rh90, B3, MNF3, B74, MNF4, rh25, TCW, B24, EVI, B7, B76, B4 ND43, rh75, B345, and rh60	20
ICESat-2	rh85, rh90, rh75, rh98, rh25, rh100, and rh60	7
Landsat 9	B7, B3, MNF2, MNF4, MNF3, TCW, pc3, B74, ND43, B1, B6, B2, B76, B24, EVI, B4, pc2, ND67, VI3, B345, TCB, and SLAVI	22
ICESat-2 and Landsat 9	rh85, rh90, rh75, B3, MNF2, B7, rh98, rh60, rh100, B76, MNF3, rh25, TCW, MNF4, B24, B74, EVI, ND43, B4, B1, ND67, pc3, and rh50	23

). By combining GEDI, ICESat-2, Landsat 9, and SRTM data, respectively, three models based on the single data source, one model based on GEDI, Landsat 9, and SRTM data source combination, and one model based on ICESat-2, Landsat 9, and SRTM data source combination were established, respectively, using the random forest regression algorithm. At last, the accuracy of the six inversion models is verified using airborne LiDAR data.

2.3.1. Importance Analysis of Feature Variables

The random forest algorithm can reduce the dimensionality of characteristic variables, thereby improving the computational efficiency of machine learning by screening variables [43]. Numerous studies have extensively shown the effectiveness of random forests in predicting forest attributes with remote sensing data [70]. In this study, the random forest method was used to analyze the importance of each characteristic variable and the first 90% of the characteristic variables with a cumulative contribution rate (the blue bars in Figure 4 represent the characteristic variables in the top 90% of the cumulative contribution rate) were screened out by calculating the importance scores of characteristic variables, and the contribution of each characteristic variable to each tree species in random forest set was evaluated. Then, the subsequent study of forest canopy height inversion was carried out.

In this study, the forest canopy height retrieval experiment was carried out based on different space-borne lidar data products after screening the feature variables by analyzing and evaluating the importance scores of the feature variables. In this study area, the GEDI model and Landsat 9 model of single data source and GEDI and Landsat 9 model of combined data source are established based on GEDI L2A products, ICESat-2 model and Landsat 9 model of single data source and ICESat-2 and Landsat 9 model of combined data source are established based on ICESat-2/ATLAS products (see Table 7). The two experiments established in this study fully evaluate the actual performance and benefit of Landsat 9 combined with different space-borne lidar data products for detailed forest canopy height retrieval so that the results of forest canopy height retrieval using multi-source remote sensing integration can be strictly analyzed and discussed.

2.3.2. Construction of the Forest Canopy Height-Retrieval Model

The random forest algorithm is capable of handling a large number of input variables, missing data, non-linear relationships, and data with non-Gaussian distributions. Random forest regression trees are a machine learning algorithm that has been successfully employed for modeling the structure and biomass of forests at regional and global scales [32,70]. In this study, we implemented a random forest regression tree ensemble method to estimate forest height.

The random forest regression model in the study [71] was trained based on the random forest regulator function of the Scikit-leam library [72] in Python. To minimize the unexpected error caused by splitting the training set and the test set, the final screened data set of the study area was split into the training set and the validation set using a ten-fold cross-validation method. The principle of cross-validation is to divide the data set into ten subsets and use one subset as the validation set and the rest as the training set. We change the verification subset in different iterations. The evaluation scores for each iteration are averaged to obtain a more robust evaluation score for the model. A 10-fold cross-validation method was used to find the optimal parameters (see

Table 8. Optimal parameters of the random forest model from different data sources.

DATA GROUP	Parameter
	n_estimators	max_depth	min_samples_split	min_samples_leaf	max_features
GEDI	288	15	1	5	0.96
Landsat 9	312	17	2	6	0.89
GEDI and Landsat 9	300	15	1	6	0.97
ICESat-2	280	12	1	6	0.92
Landsat 9	310	16	2	7	0.89
ICESat-2 and Landsat 9	200	14	1	5	0.90

) [73,74] and to build the optimal model.

According to the above information. The random forest canopy height inversion model was constructed using single and combined data sources. Table 8 lists the final parameter configurations for each model.

2.4. Accuracy Assessment

Utilizing the DTM and CHM data in NEON as reference data, terrain and canopy height were evaluated.

The correlation coefficient (R), root-mean-square error (RMSE), mean absolute error (MAE), and relative root-mean-square error (rRMSE) are calculated to measure the accuracy of each dataset in the model.

$$\mathrm{R}=\frac{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)\left({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right)}}{\sqrt{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}}\right)}^2\cdot {\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right)}^2}}}}$$

(1)

$$\mathrm{MAE}=\frac{1}{\mathrm{n}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{y}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}\right|}$$

(2)

$$\mathrm{RMSE}=\sqrt{{\frac{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{y}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}\right)}}{\mathrm{n}}}^2}$$

(3)

$$\mathrm{rRMSE}=\frac{100\%}{\overline{\mathrm{y}}}\sqrt{{\frac{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{y}}_{\mathrm{i}}-{\mathrm{x}}_{\mathrm{i}}\right)}}{\mathrm{n}}}^2}$$

(4)

where $\mathrm{n}$ represents the total number of samples, ${\mathrm{y}}_{\mathrm{i}}$ denotes the verification value, ${\mathrm{x}}_{\mathrm{i}}$ is the prediction value, $\overline{\mathrm{y}}$ represents the average of the verification value, and $\overline{\mathrm{x}}$ represents the average of the prediction value.

3. Results

3.1. Utilizing GEDI L2A Product for Forest Canopy Height Retrieval

In this study area, 5398 valid footprints screened using GEDI L2A products were used for accuracy validation to assess the ability of GEDI L2A data to accurately estimate forest canopy height. In addition to the six algorithms, GEDI L2A provides a default optimal algorithm for relative canopy height [75]. The study used a default algorithm to obtain a set of relative altitude metrics (RH90, RH92, RH94, RH96, RH98, and RH100) for each footprint in the study area and assess the correlation with the airborne radar RH90.

Table 9. The accuracy of forest canopy and reference airborne lidar forest canopy was assessed using default algorithmic methods in GEDI L2A products.

Relative Height (RH)	R	MAE/m	RMSE/m	rRMSE	N
RH90	0.74	4.37	6.16	28.21%	5398
RH92	0.75	3.93	5.67	25.94%	5398
RH94	0.77	3.55	5.24	23.98%	5398
RH96	0.78	3.27	4.96	22.71%	5398
RH98	0.78	3.22	4.94	22.59%	5398
RH100	0.78	4.00	5.63	25.76%	5398

shows the correlation results between the two. Table 9 illustrates that the retrieval accuracy of forest canopy height increases gradually from RH90 to RH98 and then decreases from RH98 to RH100. It is estimated that RH98 of GEDI L2A can be used to retrieve canopy height. This is the same result as the GEDI forest canopy height values selected from the literature [76].

We extracted values from six different sets of algorithms for RH98 data in GEDI products to investigate the correlation between data from these six sets of RH98 algorithms and airborne data for reference forest heights RH90. The results are shown in Figure 5. Among the different algorithm groups, a4 and a5 have the worst accuracy of canopy height inversion, with RMSE values as high as 9.56 m. The a2 algorithm has the highest accuracy (R = 0.78, MAE = 3.19, RMSE = 4.94 m, and rRMSE = 22.60%). Therefore, the subsequent use of an a2 algorithm data inversion model construction.

3.2. Utilizing the ICESat-2 ATL08 Product for Forest Canopy Height Retrieval

We extracted a set of relative height metrics from each footprint in the study area in ICESat-2 ATL08 data: RH75, RH85, RH90, RH98, and RH100. For ATL03 data, we obtained signal photons by matching them to corresponding ATL08 data products, extracted forest height indicators for each 25-meter strip from these signal photons, and assessed correlation with reference forest height airborne radar RH90. Figure 6 shows the accuracy of ICESat-2/ATLAS-derived RH indicators versus reference forest heights. Of all relative height indicators, RH90 performed best in estimating forest height. Therefore, RH90 derived from ICESat-2 was selected as the forest height index. Figure 7 shows a scatter plot of the RH90 index versus reference forest height for ICE-Sat-2/ATLAS day and night data. The results show that the data collected at night have a strong ability to estimate forest height. In contrast, ICESat-2 data collected during the day are unsuitable for forest height retrieval because of their low accuracy. Therefore, the inversion model is constructed by using the night acquisition data.

3.3. Utilizing the Random Forest Regression Model for Forest Canopy Height Retrieval

After preprocessing different space-borne lidar data products in the study area, we use a random forest regression algorithm to construct a forest canopy height-retrieval model for different space-borne lidar data products in the study area.

Firstly, 5398 valid laser footprints were screened from GEDI data. To optimize the efficiency and accuracy of model training, refer to Table 6 and Table 7 mentioned above for the selected model input characteristic variables and parameters. In the experimental area, there are two models based on single data sources, the GEDI model and Landsat 9 model, and there is a model based on a combination of multiple data sources, the GEDI and Landsat 9 model. The inversion results are shown in Figure 8. In models using a single data source, GEDI outperformed Landsat 9 with an increase of 0.07 in R, a decrease of 0.69 in MAE, and a decrease of 0.51 m in RMSE. This is because GEDI data provide relative heights of forest canopy structure parameters that were more important than Landsat 9 feature variables in the previously performed variable importance analysis. Among models that combine data sources, we observe significant accuracy of the GEDI, Landsat 9, and SRTM combined models. Compared to the least accurate Landsat 9 model (R = 0.71, MAE = 3.45 m, RMSE = 4.86 m, and rRMSE = 22.12%), combined model based on GEDI, Landsat 9, and SRTM data (R = 0.92, MAE = 1.91 m, RMSE = 2.78 m, and rRMSE = 12.64%) R accuracy increased by 29.58%, MAE accuracy increased by 44.64%, RMSE accuracy increased by 42.80%, and rRMSE accuracy increased by 42.86%. Compared with the optimal model in the literature [43], the MAE accuracy of the combined model based on GEDI, Landsat 9, and SRTM data in this study is improved by 13.18%, and RMSE accuracy is improved by 10.61%. Overall, accuracy has been greatly improved.

Secondly, 1690 valid laser footprints were screened from ICESat-2/ATLAS data. To optimize the efficiency and accuracy of model training, refer to Table 6 and Table 7 mentioned above for the selected model input characteristic variables and parameters. In the experimental area, there are two models based on single data sources, the ICESat-2 model and Landsat 9 model, and one model based on a combination of multiple data sources, the ICESat-2 and Landsat 9 model. The inversion results are shown in Figure 9. ICESat-2 outperformed Landsat 9 in models using a single data source with an increase of 0.24 in R, a decrease of 0.88 in MAE, and a decrease of 0.98 m in RMSE. This is because ICESat-2 data provide relative heights of forest canopy structure parameters that were more important than Landsat 9 feature variables in the previously performed variable importance analysis. Among models that combine data sources, we observe significant accuracy of ICESat-2, Landsat 9, and SRTM combined models. Compared to the least accurate Landsat 9 model (R = 0.46, MAE = 3.42 m, RMSE = 5.02 m, and rRMSE = 22.21%), combined model based on ICESat-2, Landsat 9, and SRTM data (R = 0.89, MAE = 1.84 m, RMSE = 2.54 m, and rRMSE = 10.75%) R accuracy increased by 93.48%, MAE accuracy increased by 46.20%, RMSE accuracy increased by 49.40%, and rRMSE accuracy increased by 51.60%. Compared with the optimal model in the literature [43], the MAE accuracy of the combined model based on ICESat-2, Landsat 9, and SRTM data in this study is improved by 16.36%, and the RMSE accuracy is improved by 18.33%. Overall, accuracy has been greatly improved.

4. Discussion

We analyzed and evaluated different space-borne lidar data products against reference canopy heights (90th percentile of airborne radar canopy height data). For GEDI data, we compared the reference canopy height to the data generated using the default algorithm obtained from the GEDI waveform return metrics. We analyzed six relative height (RH) measurements from RH90 to RH100 and found that RH98 had the highest correlation with reference canopy height. This finding is consistent with the findings of the study [76]. However, contrary to the results of Potapov et al. [38], RH95 was identified by Potapov et al. as the most representative measure of true forest height. A possible explanation for this discrepancy is that Potapov et al.’s study included not only forests but also other land cover types such as shrubs, grasslands, wetlands, cropland, and tundra in their canopy height estimates, resulting in an underestimation of forest height. GEDI L2A version 2 product data provides six different algorithms to improve the accuracy of GEDI L2A data by adjusting three elements of the waveform signal (smoothing width, start threshold, and end threshold) [75]. Comparing the RH98 data with reference canopy height data, it is found that Algorithm 4 and Algorithm 5 have insufficient accuracy in canopy height inversion, while Algorithm 2 has the highest accuracy. Table 5 shows that Algorithm 4 has the highest waveform start and end thresholds, while Algorithm 5 has the lowest waveform start and end thresholds. Therefore, it can be inferred that both too-short and too-long waveform lengths will result in reduced accuracy of canopy height retrieval.

For ICESat-2/ATLAS data, we extracted five forest height indicators for each 25 m band from signal photons obtained by matching ATL03 to corresponding ATL08 data products and validated them against reference canopy heights, finding the highest correlation between RH90 measurements and reference canopy heights. We then assessed the correlation between RH90 indicators and reference canopy heights in ICESat-2/ATLAS data collected during the daytime and at night in the study area. The results show that the data collected at night have a strong ability to estimate forest height. In contrast, the ICESat-2 data collected during the day are unsuitable for forest height retrieval due to their low accuracy, which may be caused by the fact that the data collected at night are not affected by excessive solar noise [77].

Multivariate analysis was performed on the characteristic variables of different space-borne lidar data products in the study area using a random forest model to determine the relative importance of each variable and select the final set of variables that fit the model. The results show that topographic factors are key variables for canopy height inversion (see Figure 10). This finding is consistent with discussions in the literature [22]. Terrain and slope are identified as the primary influencers affecting terrain height-retrieval accuracy, while canopy height emerges as the foremost error factor in canopy height retrieval. Therefore, terrain factors were added in subsequent model retrieval to improve the accuracy. The Landsat satellite series does not allow for direct measurement of forest canopy height. However, relevant feature variables can be obtained through multispectral data, which, when combined with other related variables, can estimate forest canopy height. The accuracy of the model hinges on the presence and quality of reference data, the temporal and spatial coherence of optical data, and the geographic scale of its application [38]. The literature [69] confirms the significant role of vegetation indices in estimating forest canopy height, establishing a relationship between the two using Landsat remote sensing imagery. In this study, we conducted a multivariate analysis using random forest to select the feature variables obtained from multispectral data.

Six models were established based on different space-borne lidar data products using a random forest regression algorithm to retrieve canopy height in the experimental area. In GEDI data and ICESat-2/ATLAS data, the multivariate analysis in Figure 4 shows that the importance scores for the characteristic variables of Lansat 9 data are relatively lower than the relative heights of GEDI and ICESat-2. The accuracy of canopy height retrieved using the Landsat 9 model is lower than the GEDI model and ICESat-2 model. This may be because canopy height is the main error factor in canopy height retrieval, and multispectral data cannot be measured directly.

In this study, we effectively utilized random forest regression models in machine learning to estimate forest canopy height by integrating GEDI L2A data from space-borne lidar, Landsat 9 data from optical remote sensing images and topographic factors, and ICESat-2 data from space-borne lidar, Landsat 9 data from optical remote sensing images and topographic factors, respectively. The integrated inversion results of multi-source remote sensing data are satisfactory. Compared with the results of tree height inversion using the BP neural network for multi-source remote sensing data integration in the literature [43], the inversion accuracy of two multi-source remote sensing data integration in this study is improved by 10.61% and 18.33%, respectively. The results highlight the superiority of forest canopy height inversion using the random forest regression model in this study. However, the inherent uncertainty of the random forest regression model makes it prone to overfitting at larger decision tree depths, while smaller depths may result in underfitting. Therefore, it is necessary to further strengthen parameter optimization in future research.

5. Conclusions

This study evaluated the effectiveness of multi-source remote sensing data integration and machine learning methods to improve canopy height-retrieval capabilities. In boreal temperate forest ecosystems, Landsat 9 data were combined with GEDI and ICESat-2 data representing space-borne lidar, respectively, to generate forest canopy height inversion models using random forest regression algorithms. In addition, the influence of various factors on the retrieval accuracy of canopy height is analyzed.

The study based on data from different space-borne lidar products shows that (1) the accuracy of canopy height estimation from ICESat-2 and GEDI data varies with data acquisition schemes. Among the six algorithms for GEDI L2A second edition data, the a2 algorithm estimates CHM with the best accuracy. The ICESat-2 ATL08 sixth edition data acquisition accuracy at night is better than the data collected during the day. (2) A random forest algorithm was used to rank the importance of the characteristic variables. It was found that terrain factors were important factors affecting the accuracy of canopy height retrieval, and canopy height had the greatest influence on canopy height retrieval. (3) The random forest model of multi-source remote sensing data integration is always better than the single-source model, and the Landsat 9 model is lower than the GEDI model and ICESat-2 model. The RMSE of GEDI and Landsat 9 and SRTM is 2.78 m. Compared with the GEDI model, the precision of the GEDI and Landsat 9 and SRTM models is improved by 35.95% and 42.80%, respectively. The RMSE of ICESat-2 and Landsat 9 and SRTM is 2.54 m. Compared with the ICESat-2 model, the precision of ICESat-2 and Landsat 9 and SRTM is 37.13% higher, and compared with the Landsat 9 model, the precision of ICESat-2 and Landsat 9 and SRTM is 49.40%.

In summary, this study is a preliminary exploration of the integration method of optical remote sensing Landsat 9 data with space-borne lidar products GEDI and ICE-Sat-2 of different models. The results will provide valuable insights into the practical performance and benefits of future use of Landsat 9 data in combination with other data for detailed forest canopy height retrieval, highlighting the enormous benefits of combining multispectral data and topographical factors when analyzing forest ecosystems. They provide important insights for future efforts in this field, highlighting in particular the effectiveness of combining different space-borne lidar products and Landsat 9 data separately with advanced machine learning algorithms. Nevertheless, the study acknowledges certain limitations, including limitations on site selection. Future research should focus on broader, more diverse global site selection and data accumulation, ideally incorporating GEDI and ICESat-2 data alongside other datasets like the Sentinel-3 and GaoFen-7 data. This expanded approach would likely yield more comprehensive insights for regional-scale forest monitoring.