Large-scale forest stand height mapping in the northeastern US and China using L-band spaceborne repeat-pass InSAR and GEDI lidar data

This paper presents a global-to-local fusion approach combining spaceborne synthetic aperture radar (SAR) interferometry (InSAR) and lidar to create large-scale mosaics of forest stand height. The forest height estimates are derived based on a semi-empirical InSAR scattering model, which links the forest height to repeat-pass InSAR coherence magnitudes. The sparsely yet extensively distributed lidar samples provided by the Global Ecosystem Dynamics Investigation (GEDI) mission enable the parameterization of the signal model at a finer spatial scale. The proposed global-to-local fitting strategy allows for the efficient use of lidar samples to determine the adaptive model at a regional scale, leading to improved forest height estimates by integrating InSAR–lidar under nearly concurrent acquisition conditions. This is supported by fusing the second generation of the Advanced Land Observing Satellite (ALOS-2) and GEDI data at several representative forest sites. This approach is further applied to the open-access ALOS InSAR data to evaluate its large-scale mapping capabilities. To address temporal mismatch between the GEDI and ALOS acquisitions, disturbances such as deforestation are identified by integrating ALOS-2 backscatter products and GEDI data. A modified signal model is further developed to account for natural forest growth over temperate forest regions where the intact forest landscape, along with forest height, remains quite stable and only changes slightly as trees grow. In the absence of detailed statistical data on forest growth, the modified signal model can be well approximated using the original model at the regional scale via local fitting. To validate this, two forest height mosaic maps based on the open-access ALOS-1 data were generated for the entire northeastern regions of the US and China with total area of 18 and 152 million ha, respectively. The validation of the forest height estimates demonstrates improved accuracy achieved by the proposed approach compared to the previous efforts, i.e., reducing from a 4.4 m RMSE at a few-hectare pixel size to 3.8 m RMSE at a sub-hectare pixel size. This updated fusion approach not only fills in the sparse spatial sampling of individual GEDI footprints, but also improves the accuracy of forest height estimates by 20 % compared to the interpolated GEDI maps. Extensive evaluation of forest height inversion against Land, Vegetation, and Ice Sensor (LVIS) lidar data indicates an accuracy of 3–4 m over flat areas and 4–5 m over hilly areas in the New England region, whereas the forest height estimates over northeastern China are best compared with small-footprint lidar validation data even at an accuracy of below 3.5 m and with a coefficient of determination (R²) mostly above 0.6. Given the achieved accuracy for forest height estimates, this fusion prototype offers a cost-effective solution for public users to obtain wall-to-wall forest height maps at a large scale using freely accessible spaceborne repeat-pass L-band InSAR (e.g., forthcoming NISAR) and spaceborne lidar (e.g., GEDI) data. These products are available via https://doi.org/10.5281/zenodo.11640299 (Yu and Lei, 2024).

1 Introduction

Forests play a crucial role in the terrestrial ecosystem as they serve as one of the largest terrestrial carbon pools (Pan et al., 2011). As identified by IPCC v6 (Masson-Delmotte et al., 2021) and international forest monitoring efforts such as the United Nation's REDD+ program (Angelsen, 2009), large-scale (e.g., state, continental, and global) forest height products are desired to quantify carbon storage in forested resources due to their close relationship to aboveground biomass (AGB). These products also help to determine forests' roles in climate change mitigation and biodiversity conservation (Houghton et al., 2009). In this work, “forest height” or “forest stand height” (FSH) is referred to as the medium-footprint (25 m) lidar-determined relative height at the 98th percentile (RH98) as measured by NASA's GEDI lidar mission on board the International Space Station (ISS).

Satellite-based remote sensing represents a cost-effective method to investigate biophysical parameters of forests. Commonly used remote sensing methods include optical and microwave imaging observations, such as passive optical sensors including Landsat series (Loveland and Dwyer, 2012); lidar missions including ICESat-1/2 (Schutz et al., 2005; Abdalati et al., 2010) and Global Ecosystem Dynamics Investigation (GEDI; Dubayah et al., 2020) missions; and synthetic aperture radar (SAR) systems such as Japanese Aerospace Exploration Agency (JAXA)'s ALOS-1/2 (Rosenqvist et al., 2007, 2014), Sentinel-1 (Torres et al., 2012), and TanDEM-X (Krieger et al., 2007). Lidar and SAR are promising for capturing the internal vertical structure of forests: lidar is fundamentally sensitive to structural details, while radar detects the three-dimensional distribution of vegetation elements (Ulaby et al., 1990). The backscatter information from a single SAR image can be used for inferring the AGB (Santoro et al., 2021) despite the fact that the actual vertical information remains undetermined. As an extension of SAR backscatter observation, SAR interferometry (InSAR) provides direct information related to the vertical forest structure (Treuhaft and Siqueira, 2000). Spaceborne InSAR can operate by either single-pass bistatic interferometer (e.g., TanDEM-X missions; Krieger et al., 2007) or repeat-pass InSAR (e.g., ALOS-1/2 L-band and Sentinel-1 C-band missions). The short wavelength operated in former satellites may restrict its sensing capabilities over dense forests, while temporal decorrelation affects repeat-pass InSAR performance (Zebker and Villasenor, 1992; Monti-Guarnieri et al., 2020; Lavalle et al., 2012; Ahmed et al., 2011). Lidar has been widely used for characterizing the forest vertical structure at a regional scale and can serve as a benchmark for calibrating inversion models and forest height estimates (Choi et al., 2023; Askne et al., 2013). Spaceborne lidar (Schutz et al., 2005) has been further developed for global ecosystem monitoring. Because of observational constraints, these measurements have been acquired based on such a spatial sampling technique that collect sparse yet extensive measurements.

For instance, NASA's GEDI mission is the first spaceborne lidar instrument designed to study ecosystems. Since 2019, GEDI has provided extensively distributed lidar waveform measurements covering nearly all global forests. These waveform observations allow for the extraction of various biophysical parameters, such as the canopy height and leaf area index. However, GEDI collects only discrete footprint measurements, spaced approximately 60 m apart in the along-track direction and 600 m apart in the cross-track direction. To overcome this limitation and extend GEDI's measurements into continuous datasets, several fusion studies have been conducted. Notable examples include efforts that incorporate radiometric information from optical sensors, such as NASA's Landsat (Potapov et al., 2021) and ESA's Sentinel-2 (Lang et al., 2022), as well as from SAR backscatter signals (Shendryk, 2022). However, relying solely on radiometric information to expand lidar observations has proven suboptimal, particularly in high-biomass regions where signal saturation occurs (Kalacska et al., 2007; Imhoff, 1995; Minh et al., 2014).

In contrast, because of its fundamental sensitivity to height and/or variations in height, the fusion of SAR interferometry and GEDI has gained much interests. For example, making joint use of TanDEM-X and GEDI data has been assessed and demonstrated for achieving wall-to-wall forest height and AGB mapping (Qi and Dubayah, 2016; Choi et al., 2023; Guliaev et al., 2021; Qi et al., 2019). Without temporal decorrelation effects, TanDEM-X data offer opportunities to leverage very high resolution observations for addressing spatially heterogeneous landscapes. However, the forest height was inverted in these studies based on the Random Volume Over Ground (RVoG) model (Treuhaft and Siqueira, 2000; Cloude and Papathanassiou, 1998), and an external constraint was induced by GEDI waveform information. That is, only the mean waveform information across the scene in these studies was used for model-based inversion, implying an underlying assumption that the forest objects over the scene share a similar vertical structure. This may lead to a degraded performance when dealing with spatially heterogeneous forests. To address this, Qi et al. (2025) proposed a regional postprocessing correction model to refine suboptimal height estimates, while Hu et al. (2024) exploited local ICESat-2 lidar information, using regional polynomials and an adaptive window, to estimate equivalent forest phase centers under homogeneous forest and terrain conditions. Additionally, a potential limitation of TanDEM-X observations is the insufficient penetration capability over dense forests due to the short wavelength of the X-band (∼ 3.1 cm) (Kugler et al., 2014). This underlines the need for longer-wavelength SAR systems (e.g., L-band) to enhance sensitivity in dense forest environments.

Temporal decorrelation has been a widely studied topic in InSAR research (Rocca, 2007; Ahmed et al., 2011; Bhogapurapu et al., 2024). Zebker and Villasenor (1992) proposed a Gaussian model to analyze oceanic scenarios, while Monti-Guarnieri et al. (2020) summarized the signal models tailored for vegetated scenarios. Askne et al. (1997) introduced a coordinate dependence of the vertical motion profile to analyze InSAR temporal decorrelation effects caused by wind. Building upon the well-known RVoG model, several signal models have been developed to explicitly incorporate temporal decorrelation effects (Lavalle et al., 2012; Papathanassiou and Cloude, 2003; Lei et al., 2017a).

This study employs the RVoG-based temporal decorrelation model (Lei et al., 2017a) to invert forest height. The model-based inversion was first demonstrated using a relatively small airborne lidar reference strip to generate a forest height mosaic over a two-state region in the northeastern US (Lei et al., 2018). Given the limited availability of lidar datasets at that time, scene-wide constant model parameters for the relationship between the InSAR temporal decorrelation and lidar observations were assumed, and the overlapping area between InSAR scenes was used to propagate the lidar information throughout the adjacent InSAR scenes. The advent of NASA's Global Ecosystem Dynamics Investigation (GEDI) mission has since enhanced this methodology by integrating local GEDI samples directly into the inversion framework (Lei and Siqueira, 2022; Yu et al., 2023). This integration enables spatially adaptive calibration of model parameters, overcoming prior limitations of constant parameter assumptions and improving inversion accuracy in heterogeneous forest land covers.

This paper further removes the assumption of a spatially constant temporal change model that was made in the previous efforts and develops a new inversion approach based on a two-stage (global-to-local) inversion strategy. By efficiently leveraging regional GEDI samples, this approach calibrates a semi-empirical, semi-physical repeat-pass InSAR model at a finer spatial scale, substantially improving forest height inversion accuracy. The method assumes that the temporal decorrelation model remains spatially invariant at the regional scale while permitting variability in forest height observations within those regions. This approach is validated by fusing ALOS-2 InSAR and GEDI data acquired under nearly concurrent conditions. Furthermore, the approach is applied to the open-access ALOS InSAR data for evaluating its large-scale mapping capability. To address the temporal mismatch between the ALOS and GEDI acquisition, forest disturbance can be detected by fusing SAR backscatter and lidar data under nearly concurrent conditions. Furthermore, a modified model is developed to account for the natural growth of forests over temperate forest regions where the intact forest landscape (IFL) exhibits slow changes in height. Without available forest growth data, the modified signal model can be well approximated using the original signal model at the regional scale through local fitting. Two 30 m gridded forest height mosaics were generated for the northeastern regions of US and China. Validation of the generated forest height mosaics against extensive airborne lidar observations demonstrates enhanced inversion accuracy at sub-hectare pixel size. The key contribution of this paper lies in the use local GEDI information for radar–lidar data fusion, enabling large-scale and efficient forest height mapping using open-access spaceborne data, such as GEDI and forthcoming NISAR (Siqueira et al., 2024; Kellogg et al., 2020) data.

This paper is structured as follows: Sect. 2 describes the study areas and remote sensing datasets, Sect. 3 details the proposed inversion framework, and Sect. 4 validates forest height estimates across diverse forest sites in the northeastern US and China. Section 5 discusses the implications and limitations of the methodology, and Sect. 8 concludes this study.

2 Study area and datasets

2.1 Study area

This paper focuses on the northeastern regions of the US and China. These regions contain transitional forests composed of both coniferous and broad-leaved species. As shown in Fig. 1, the New England region in the northeastern US (including the states of Maine, New Hampshire, Vermont, Massachusetts, Connecticut, and Rhode Island) is selected for the generation and validation of the large-scale mosaic of forest height (covering a total area of 18 million ha) due to the availability of ample airborne lidar datasets. Forests in this area are primarily dominated by coniferous forests (red pine, balsam fir, hemlock etc.) and northern hardwoods (maple, oak, beech, etc.). These forests exhibit stable intact forest landscapes (IFLs) and canopy heights, with changes driven primarily by natural growth rather than anthropogenic disturbance (Riofrío et al., 2023).

Figure 1 Study area and validation sites for the New England region in the US. The generated forest height mosaic map covers the states of Maine, New Hampshire, Massachusetts, Vermont, Connecticut, and Rhode Island. The inversion results are validated against the large-footprint (25 m) LVIS data acquired in either 2009 or 2021 over the validation sites (denoted by red dot markers). At the White Mountain National Forest (WMNF) site, small-footprint GRANIT lidar data are also used for validation after reprocessing into equivalent RH98 metric maps. The features of the validation sites are summarized in Table 1. The optical basemap is from © Microsoft Bing Maps.

Another large-scale forest height mosaic is also generated over the northeast of China with a total area of 152 million ha. As shown in Fig. 2, the forest height mosaic for China covers five provinces: Hebei, Jilin, Liaoning, Inner Mongolia, and Heilongjiang. Note that Jilin and Inner Mongolia provinces were not fully covered as only forested areas within the GEDI observation coverage (< 51.6° N) are addressed. The forests in northeastern China can be primarily grouped into four primary regions: (1) deciduous coniferous forest region located at the northernmost parts of Inner Mongolia and Heilongjiang provinces, (2) temperate mixed-forest region (comprising evergreen coniferous and deciduous broad-leaved species) primarily distributed in Heilongjiang and Jilin provinces, (3) northern temperate mixed-forest subregion situated in Liaoning province, and (4) temperate steppe region located partly in Hebei province and partly in Inner Mongolia province.

Figure 2 Study area and validation sites in northeastern China. Five provinces are covered in the generated forest height mosaic: Jilin, Liaoning, Hebei, Heilongjiang, and Inner Mongolia. The performance of the forest height inversion is assessed by comparison with small-footprint (0.5–1 m) lidar data at the validation forest sites (indicated by the red diamond markers) in each province. The features of the validation sites are summarized in Table 2. The optical basemap is from © Microsoft Bing Maps.

The northeastern regions of the US and China were selected as study areas for two key reasons. First, both regions provide access to extensive airborne lidar datasets: NASA's Land, Vegetation, and Ice Sensor (LVIS) (Blair et al., 1999) in the US and small-footprint airborne lidar data in China. Second, the performance of forest height inversion can be assessed in a unique way: the New England region offers abundant GEDI calibration sites, while the northeastern part of China is situated at a comparable region at a similar latitude but without dedicated GEDI calibration sites.

Several experimental and validation sites are selected across the northeastern regions of the US and China, with their relevant information summarized in Tables 1 and 2, respectively. In the New England region, validation is conducted at various sites including the Howland Forest site in Maine, the Harvard Forest site in Massachusetts, the White Mountain National Forest (WMNF) site in New Hampshire, the Green Mountain National Forest (GMNF) site in Vermont, and the Naugatuck State Forest (NSF) site in Connecticut. These forest validation sites are covered by medium-footprint (25 m) LVIS data acquired in either 2009 or 2021. Particularly, the forest height inversion over the WMNF site was evaluated using GRANIT airborne laser scanning (ALS) data acquired in 2011 (Haans et al., 2009), with the canopy height product extracted from the waveform data at a raster sampling spacing of 2 m.

Table 1 The forest validation sites covered by the airborne lidar observation in New England, US.

Regarding northeastern China, evaluation was performed at one forest site in each province: the Mengjiagang Forest site in Heilongjiang province, the Dagujia Forest site in Liaoning province, the Saihanba Forest site in Hebei province, the Hubao National Park in Jilin province, and the Genhe Forest Bureau in Inner Mongolia province. Validation of the forest height product across all the forest validation sites in China was done by comparisons with the small-footprint (0.5–1 m) ALS data.

To provide a preliminary assessment of forest disturbances in these two regions, forest disturbance maps derived from global forest change products (Hansen et al., 2013) are shown in Figs. 3 and 4.

Table 2 The forest validation sites covered by the ALS validation data in northeastern China.

Figure 3 Forest disturbance map of the New England region (2007–2023) derived from the Global Forest Change dataset (Hansen et al., 2013). The binary classification distinguishes undisturbed areas (0) from disturbed areas (1) within the period. The optical basemap is from © Microsoft Bing Maps.

Figure 4 Forest disturbance map of northeastern China (2007–2023) derived from the Global Forest Change dataset (Hansen et al., 2013). The binary classification distinguishes undisturbed areas (0) from disturbed areas (1) within the period. Province names are shown in green on the central map. The optical basemap is from © Microsoft Bing Maps.

2.2 Spaceborne and airborne remote sensing datasets

The freely accessible L-band InSAR data from the Japanese Aerospace Exploration Agency (JAXA)'s Advanced Land Observing Satellite-1 (ALOS-1) mission were used for generating the InSAR correlation observations. In addition, some radar data from ALOS-2 (a follow-up mission of ALOS-1) were also employed as case studies for validating the fusion of radar–lidar under concurrent conditions. Furthermore, the spaceborne lidar waveform-based metrics (RH98 as a proxy for forest height) from the GEDI mission were used for parameterizing the temporal decorrelation model.

2.2.1 Spaceborne InSAR datasets

Global fine beam dual-polarization (FBD) SAR images with a spatial grid of 10 × 3 m (for range and azimuth directions) were collected by the ALOS-1 satellite from 2007 to 2010, with a repeat cycle of 46 d. To generate a forest height mosaic, 100 cross-polarized ALOS-1 InSAR scenes (identified by one pair of frame and orbit numbers) were processed to cover the New England region in the US, whereas more than 600 InSAR scenes were processed for covering the northeast of China. The InSAR preprocessing was done by the Jet Propulsion Laboratory (JPL)'s InSAR Scientific Computing Software (ISCE) software (Rosen et al., 2012). It was reported by Lei and Siqueira (2014) that the ALOS-1 InSAR observations acquired during the summer/fall time frame of 2007 and 2010 tended to have higher InSAR coherence. In practical processing, multiple cross-polarized interferograms during the lifetime of ALOS-1 (2007–2011) were formed and processed for each scene based on different combinations of acquisition dates, allowing for the identification of the best InSAR pair for each ALOS-1 observation scene.

As a follow-on mission to ALOS-1, ALOS-2 had a shorter revisiting period (14 d), resulting in a better InSAR correlation behavior. The acquisition started from 2016, allowing for nearly concurrent observations with respect to GEDI samples. However, the acquisition strategy and limited access to the high-resolution dual-polarized strip-map data have made it more difficult to form proper InSAR pairs and perform large-scale mapping. The grid of ALOS-2 images in FBD mode is at a grid of 8 × 4 m for range and azimuth directions. A window size of 4 × 8 looks is used for the coherence estimation, resulting in an average pixel size of 30 m.

In this study, the use of limited ALOS-2 data is devoted to demonstrating the proposed approach in the ideal case. The large-scale mapping capabilities were demonstrated using free-access ALOS-1 data over the temperate forest regions. It is also noted that there is a time discrepancy between the acquisition dates of ALOS-1 and GEDI. This discrepancy is addressed using the twofold solution as detailed in Sect. 3.3. As for the abrupt discrepancy due to forest disturbance (e.g., logging, deforestation, and fire) that usually results in no/short vegetation with small backscatter values, replacing the InSAR inverted forest estimates with those derived from the appropriate ALOS-2 backscatter mosaic map for short vegetation (as shown in Sect. 3.1) can detect the disturbed forest areas, so the study area is mainly concentrated on temperate and boreal forests, the heights of mature temperate forests (intact forest landscape) remain almost stable with slight changes. Nevertheless, a simulation in Sect. 3.4.3 shows that the approach of this study can approximate the height of the forests subject to natural growth at a regional scale. In other words, all the InSAR-based height inversions are calibrated to the acquisition time of GEDI and thus best compared with the concurrent airborne lidar validation data.

2.2.2 Spaceborne lidar datasets

As the first spaceborne lidar mission to characterize ecosystem structure and its dynamics, the NASA's GEDI mission was launched in December 2018. GEDI provides near-global measurements of forest structure metrics from 51.6° S to 51.6° N until 2024. With three lasers mounted, eight parallel tracks of samples at a footprint of 25 m are simultaneously collected. The spatial separation between samples during one data take is 60 m in the along-track direction and 600 m in the cross-track direction. The GEDI RH98 metric is selected as an appropriate proxy for indicating forest canopy height within each footprint because it has a lower sensitivity to errors as compared to the RH100 metric (Hofton et al., 2020). After filtering out GEDI samples with lower penetration sensitivities (e.g., 95 % sensitivity, 50 m maximum elevation difference between GEDI and TanDEM-X measurements), the remaining L2A version 2 GEDI samples are used for calibrating the inversion model.

2.2.3 Airborne lidar data

Medium-footprint lidar data

A significant amount of airborne full-waveform lidar data was collected across the US using the LVIS sensor. The lidar data were processed into RH98 maps with a 25 m grid. Specifically, lidar data over the Howland Forest site in Maine were obtained in 2009. LVIS data acquired in 2021 for GEDI calibration cover all the other forest sites in the New England region, which were classified into four parts based on the state boundaries of New Hampshire, Vermont, Massachusetts, and Connecticut.

Small-footprint lidar data

In some sites, small-footprint lidar data have to be used for validation. The validation at the WNMF site utilized GRANIT lidar data, which has a 2 m footprint and was acquired in 2011. All the forest heights in northeastern China were validated using small-footprint lidar data. The airborne lidar data, with an average point density of 6 pts m⁻² over Hubao National Park, were acquired using an airborne lidar system owned by the Chinese Academy of Forest Inventory and Planning. The observations covering all other forest sites in the northeastern China were acquired during 2017–2022 using the airborne remote sensing system developed by the Chinese Academy of Forestry (Pang et al., 2016), which has an average point density of 12 pts m⁻². It should be noted that validating the forest height estimates against airborne LVIS observations is not straightforward as the footprint of these airborne data is much smaller than the footprint of GEDI. To address this footprint difference, an equivalent RH98 metric (referred to as ERH98 hereafter) needs to be extracted at the position of the 98th percentile of the lidar waveform or from the histogram formed by high-resolution canopy height model (CHM) estimates within a matching footprint size. Following this procedure, all small-footprint lidar data were reprocessed to generate forest height estimates based on the ERH98 metric using a 25 m footprint

Figure 5 Block diagram of the workflow for generating forest height mosaics.

Figure 6 An illustrative example of the processing steps at the Howland Forest site: (a) the input ALOS-1 coherence magnitude map, (b) the GEDI RH98 samples, (c) the forest height estimates based on InSAR coherence information, (d) the backscatter-based height estimates, (e) the final forest height map after replacing the estimates of short trees in (c) with the collocated pixels in (d). Panel (f) shows the airborne LVIS lidar validation data.

2.2.4 Forest and non-forest maps

Non-forest areas including waterbodies and urban areas were masked out using the 2021 National Land Cover Database (Homer et al., 2015) and the 2021 ESRI Global Land Cover Map (Karra et al., 2021).

2.2.5 Backscatter mosaic map

This study used the global radar backscatter products generated by JAXA using ALOS-1/2 FBD images (Shimada et al., 2016) after radiometric and geometric calibration (including slope effects correction). Specifically, the global cross-polarized backscatter products from 2019 and 2020 over the northeast of the US and China were utilized to obtain height estimates of short vegetation. These 2-year products were used to account for missing data gaps and backscatter calibration inconsistencies and to best match the acquisition time of the validation airborne lidar data.

3 Methodology

This section first outlines the processing workflow and implementation steps and then provides a detailed analysis of the methodological principles.

3.1 Processing workflow

The entire processing workflow is illustrated in Fig. 5. Figure 6 provides a zoomed-in example of global-to-local forest height inversion over the Howland Forest site. The standard interferometric preprocessing (including coregistration, topographic-phase compensation, interferometric formation, and geocoding steps) for multiple ALOS-1 InSAR pairs is performed using JPL's ISCE software, in which the parallel computing capabilities of graphic processing units (GPUs) is utilized to enhance the processing efficiency, generating the InSAR coherence map as shown in Fig. 6a. After matching the GEDI samples onto the grid of InSAR scene (the collocated GEDI RH98 samples are shown in Fig. 6b), the global-to-local inversion for each scene is carried out as follows:

1. Global fitting. As illustrated in Sect. 3.3, a global fitting is first performed to obtain an initial guess of temporal parameters (S₀, C₀) (see Eq. 3) using all the available GEDI RH98 samples and corresponding InSAR observations over each InSAR scene. We noted earlier that a moderate forest disturbance may result in an overestimation of forest height for global fitting (Lei et al., 2017b), and a re-weighted iterative global fitting is instead used to remove the gross errors induced by forest disturbances after the first iteration.

2. Local fitting. A local fitting is then conducted around each GEDI sample within a local boxcar window with a constraint of a smaller searching range in the vicinity of the initial estimates, i.e., $\left(S\in \left[S_{\mathrm{0}}-r_S,S_{\mathrm{0}}+r_S\right],C\in \left[C_{\mathrm{0}}-r_C,C_{\mathrm{0}}+r_C\right]\right)$, where r_S and r_C represent the searching region for two parameters, respectively. Spatial distance-based weights (e.g., Fig. 7) are preferred in local fitting for preserving detailed information and ensuring the robust inversion over large-scale application. An expensive computational burden is implied for the local fitting so that a GPU-based implementation (Yu et al., 2019) is developed for enhancing processing efficiency.

3. Interpolation and inversion. The irregularly distributed temporal decorrelation parameters, i.e., $\left\{S\left(i,j\right),C\left(i,j\right)\right\}$, with i,j denoting the latitude and longitude for each GEDI sample, are interpolated onto the regular grid of the InSAR coherence magnitude map based on the Delaunay-triangulation-based natural neighborhood interpolation. The forest height is then inverted on a pixel-by-pixel basis using the InSAR coherence observation and the inversion model (generating wall-to-wall forest height mapping as shown in Fig. 6c).

4. Backscatter-based estimates. The forest height estimates of short vegetation in the previous inversion are replaced by the corresponding backscatter-based estimates (see Fig. 6d): the GEDI RH98 samples over bare ground and short vegetation land covers (i.e., $\mathrm{0}\le \mathrm{RH}\mathrm{98}\le \mathrm{10}$ m), along with the corresponding backscatter information, are fitted using an exponential model, which is then used to obtain backscatter-to-height estimates (Lei et al., 2019). Short trees are identified using a criterion where backscatter-derived forest height estimates fall below 10 m based on the maximum height of shrubs (Lawrence, 2013; Allaby, 2012) and empirical studies (Lei et al., 2019).

5. Mosaicking. A mosaicking approach is finally carried out to pick up the best InSAR pairs based on the pre-inversion metric and mosaicking the overlapping area using the pixels with a better goodness of fit as discussed in Sect. 3.5.

3.2 Sinc inversion model

After standard InSAR preprocessing (including coregistration, topographic-phase compensation, etc.), a complex InSAR correlation observation between two SAR images can be derived by

$$\begin{array}{}\text{(1)}& \mathit{\gamma }={\displaystyle \frac{\langle I_{\mathrm{1}}\cdot I_{\mathrm{2}}^{*}{\rangle }_L}{\sqrt{\langle I_i\cdot I_{\mathrm{1}}^{*}{\rangle }_L\langle I_{\mathrm{2}}\cdot I_{\mathrm{2}}^{*}{\rangle }_L}}},\end{array}$$

where I₁ and I₂ represent two SAR images from the first and second acquisition, respectively, and the operator 〈⋅〉_L is used for spatial averaging of L looks. As a metric to measure similarity, the complex InSAR correlation accounts for various forms of decorrelation (Zebker and Villasenor, 1992; Gatelli et al., 1994; Treuhaft and Siqueira, 2000). If we consider a pair of SAR images with short spatial separation over forested area, InSAR coherence can be expressed as

$$\begin{array}{}\text{(2)}& \left|\mathit{\gamma }\right|=\left|{\mathit{\gamma }}_{\mathrm{SNR}}\right|\cdot \left|{\mathit{\gamma }}_{\mathrm{geo}}\right|\cdot \left|{\mathit{\gamma }}_{\mathrm{v}\mathrm{\&}\mathrm{t}}\right|,\end{array}$$

where γ_SNR represents the decorrelation induced by the radar-return signal-to-noise ratio and γ_geo denotes the geometric decorrelation due to the difference between two looking angles. After accounting for signal-to-noise ratio (SNR) decorrelation and performing common band filtering (Gatelli et al., 1994), the remaining component of correlation, γ_v&t, is only related to temporal changes and the distribution of scatterers in the vertical direction. The spatial separation for one InSAR pair is usually indicated by the interferometric vertical wavenumber κ_z (unit rad m⁻¹) (Bamler and Hartl, 1998). A small value of κ_z (less than 0.05 rad m⁻¹) is suggested by the analysis in Lei and Siqueira (2014).

For moderate and large temporal baselines, moisture-induced dielectric fluctuations and wind-induced random motion are identified as the primary factors influencing temporal decorrelation (Lavalle et al., 2012; Askne et al., 1997). Lei et al. (2017b) introduced a modified Random Volume over Ground (RVoG) model that accounts for the coupled effects of dielectric changes and random motion in spaceborne ALOS-1/2 observations. The model is further simplified by assuming negligible ground scattering under cross-polarized (HV) observations and short spatial baselines and establishes a semi-empirical formula linking the cross-polarized InSAR coherence magnitude (${\mathit{\gamma }}_{\mathrm{v}\mathrm{\&}\mathrm{t}}^{\mathrm{HV}}$) to forest height (h_v):

$$\begin{array}{}\text{(3)}& \left|{\mathit{\gamma }}_{\mathrm{v}\mathrm{\&}\mathrm{t}}^{\mathrm{HV}}\right|=S\cdot \mathrm{sinc}\left({\displaystyle \frac{h_{\mathrm{v}}}{C}}\right),\text{with}\mathrm{0}\le h_{\mathrm{v}}\le \mathit{\pi }C,\end{array}$$

where S (ranging from 0 to 1; unitless) is a parameter primarily connected to moisture-induced dielectric changes in the target and C (a positive value that has units of meters) relates to the wind-induced random motion. In practice, ground scattering remains present in cross-polarized signals, particularly for low-frequency SAR observations. However, it is typically minimal, with a ground-to-volume ratio generally lower than 0.1. This residual ground scattering can bias forest height estimates, particularly for short or very tall forest stands. It is also important to note that this formula is valid when the interferometric vertical wavenumber is below 0.15 rad m⁻¹ and is most reliable when it is smaller than 0.05 rad m⁻¹ (Lei and Siqueira, 2014), which is consistent with the acquisition geometry of ALOS-1/2 InSAR acquisition. Otherwise, the presence of volumetric decorrelation will result in a compromised inversion performance.

Figure 7 An illustrated example of distance-based weight setting for a 0.96 km (32 by 32 pixels) wide moving boxcar window.

3.3 Fusing concurrent radar–lidar acquisitions by a global-to-local inversion approach

Determination of the above temporal decorrelation model (3) has to resort to ancillary forest height data, e.g., field inventory data and lidar measurements. A global-to-local procedure is used for determination since regional forests are usually relatively homogeneous, leading to the construction of ill-conditioned observations for model determination. A global initial guess is firstly needed for constraining the solution space for subsequent fitting.

For one InSAR scene, a constant scene-wide behavior of the temporal change parameters is firstly assumed, e.g., (S_scene,C_scene). For one candidate pair of these parameters, forest height estimates $\widehat{h_{\mathrm{v}}}$ are derived by solving model (3) based on InSAR coherence magnitude. A covariance matrix between $\widehat{h_{\mathrm{v}}}$ and the ancillary forest height h_a from lidar and its eigen-decomposition are expressed by

$$\begin{array}{}\text{(4)}& \mathrm{\Sigma }=\left[\begin{array}{cc}\mathrm{cov}\left(\widehat{h_{\mathrm{v}}},\widehat{h_{\mathrm{v}}}\right)& \mathrm{cov}\left(\widehat{h_{\mathrm{v}}},h_a\right)\\ \mathrm{cov}\left(\widehat{h_{\mathrm{v}}},h_a\right)& \mathrm{cov}\left(h_a,h_a\right)\end{array}\right]=\mathbf{Q}\mathrm{\Lambda }{Q}^{-\mathrm{1}},\end{array}$$

where the 2×2 matrix Q contains the eigen vectors, the 2×2 diagonal matrix Λ comprises the eigen values. The function cov(A,B) denotes the covariance between the two vectors A and B:

$$\begin{array}{}\text{(5)}& \mathrm{cov}\left(A,B\right)={\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}{\sum }_{i=\mathrm{1}}^N\left(A-{\mathit{\mu }}_A\right)\cdot \left(B-{\mathit{\mu }}_B\right),\end{array}$$

with μ_A and μ_B being the mean values of two input vectors. Two fitting metrics, i.e., slope k and bias b, can be extracted out of the constructed covariance matrix:

$$\begin{array}{}\text{(6)}& {\displaystyle }k={\displaystyle \frac{Q_{\mathrm{21}}}{Q_{\mathrm{11}}}},\text{(7)}& {\displaystyle }b=\mathrm{2}{\displaystyle \frac{{\mathit{\mu }}_{\widehat{h_{\mathrm{v}}}}-{\mathit{\mu }}_{{\mathrm{h}}_a}}{{\mathit{\mu }}_{\widehat{h_{\mathrm{v}}}}+{\mathit{\mu }}_{{\mathrm{h}}_a}}}.\end{array}$$

By defining a figure of merit, i.e., approaching zero bias and unity slope (i.e., the 1:1 line), a pair of temporal parameters is determined by minimizing the objective function:

$$\begin{array}{}\text{(8)}& \left\{S_{\mathrm{scene}},C_{\mathrm{scene}}\right\}=\mathrm{argmin}\left(b^{\mathrm{2}}+{\left(k-\mathrm{1}\right)}^{\mathrm{2}}\right).\end{array}$$

Resolving model parameters (S_scene,C_scene) is referred to as global fitting strategy. The effectiveness of this approach was demonstrated over the northeastern US (Lei et al., 2018) with well-constructed observation vectors including sufficient samples ranging from low trees to tall trees.¹

Because temporal change factors tend to be spatially varying with respect to vegetation types and weather/climate conditions, inversion performance for solely global fitting as described above is expected to deteriorate when applied to large-scale inversion. With the availability of GEDI data, a substantial improvement is expected by expressing the decorrelation model at a finer spatial scale. A preliminary effort was made in Lei and Siqueira (2022) where the determination of (S,C) is carried out on each GEDI lidar sample h_gedi(i,j) to obtain S(i,j) and C(i,j), where i and j represent latitude and longitude for each GEDI sample. As only one sample cannot provide sufficient observations to resolve two unknowns, an external constraint on the ratio of $S\left(i,j\right)/C\left(i,j\right)$ was introduced.

Each individual GEDI measurement is subject to errors caused by artifacts such as terrain slopes (Wang et al., 2019) and systematic geolocation inaccuracies ranging from several meters to tens of meters (Tang et al., 2023). Additionally, the penetration of the lidar signal is sometimes limited within GEDI's coverage beams. InSAR coherence estimates also experience measurement uncertainty (Rodriguez and Martin, 1992; Touzi et al., 1999). To address these challenges, spatial averaging is commonly used to reduce errors for both coherence and GEDI measurements.

In this context, a local fitting step is introduced to acquire temporal change factors at a finer spatial scale using extensive GEDI information: a circular window (32 pixels, 960 m wide) is set for each GEDI sample to collect regional samples and their corresponding InSAR coherence magnitude observations. A Euclidean norm-based fitting is used since a local window tends to encompass homogeneous vegetation with similar heights so that the global k−b fitting metric is no longer robust. Assuming consistent temporal change factors within each local window, the factors can be determined by

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& \left\{S\left(i,j\right),C\left(i,j\right)\right\}=\mathrm{argmin}\sum _{(r,c)\mathit{ϵ}W}\\ & {\left(\widehat{h_{\mathrm{v}}}\left(i+r,j+c\right)-h_{\mathrm{gedi}}\left(i+r,j+c\right)\right)}^{\mathrm{2}},\end{array}\end{array}$$

where W is the local searching window in the geographic coordinate system, with r and c being local indices along the latitude and longitude direction within the local window. To preserve detailed information, weighting factors based on the distance or similarity metric are used in local fitting, as in

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& \left\{S\left(i,j\right),C\left(i,j\right)\right\}=\mathrm{argmin}\\ & {\displaystyle \frac{\sum _{\left(r,c\right)\mathit{ϵ}W}{\left[w(r,c)\left(\widehat{h_{\mathrm{v}}}\left(i+r,j+c\right)-h_{\mathrm{gedi}}\left(i+r,j+c\right)\right)\right]}^{\mathrm{2}}}{\sum _{\left(r,c\right)\mathit{ϵ}W}{\left(w\left(r,c\right)\right)}^{\mathrm{2}}}},\end{array}\end{array}$$

where the weights can be set based on spatial distances between neighboring pixels and the central pixel within the local window or adapted depending on land cover type (Deledalle et al., 2014). For computational efficiency and versatility, the spatial-distance-based weight setting, as exemplified in Fig. 7, is used for inversion in this work. The optimal local window size is determined by minimizing the root mean square error (RMSE) of forest height estimates for moderate-to-tall trees (> 10 m) across diverse window size configurations validated against independent lidar datasets. The selection of window size is always a compromise between smooth and detailed information. This window size is selected here to include enough samples for model fitting while maintaining local detailed information. As the parameters were determined on the grid of GEDI samples, a post-interpolation is needed to obtain gridded temporal parameters for matching the InSAR scene.

It should be noted that this physical scattering model (3) is established for forested areas, meaning it does not adequately address other land cover types, such as waterbodies and urban areas. Furthermore, actively managed regions (e.g., cropland) are subject to additional decorrelation effects induced by human activities. It was suggested to use L-band backscatter-based height estimates to improve the inversion in these instances. The backscatter-based inversion are derived based on an exponential model (Yu and Saatchi, 2016; Lucas et al., 2006) which can be determined using the similar global-to-local routine as the proposed approach. In practice, only the global fitting strategy is used for model parameterization as it provides accurate estimates as compared to global-to-local fitting (slightly worse), while it maintains an affordable computational expense. After that, the estimates of active land management areas in coherence-based inversion are replaced by the corresponding backscatter-based estimates.

Figure 8 Example of fusing the ALOS-2 backscatter products and GEDI to detect forest disturbances as defined by Global Forest Change products (Hansen et al., 2013) over a representative region in the New England region. Yellow pixels represent forest disturbance areas detected by backscatter-based estimates, while red pixels indicate disturbed areas as not identified by these estimates. The optical basemap is from © Microsoft Bing Maps.

Figure 9 (a) Scatterplot of forest height growth rate versus forest height derived from a comparison of LVIS forest height estimates between 2009 and 2021 at the Harvard Forest site, Massachusetts. The red line represents the fitted growth rate function of forest height in 2009 over the high-density region. (b) The fitted forest height growth rate functions at typical forest sites of other states in the New England region.

Figure 10 Simulation results comparing the original sinc model and its modified version incorporating natural forest growth at the Harvard Forest site (with forest growth rate function of $G=-\mathrm{0.0134}\cdot h_{\mathrm{v}}\left(t_{\mathrm{1}}\right)+\mathrm{0.464})$. The blue line represents the original sinc model with established temporal parameters, while the red line shows the modified model accounting for forest height growth. The dashed orange line highlights the region that cannot be characterized by the modified model. The y axis represents the coherence magnitude estimated from Eqs. (3) and (15).

3.4 Fusing non-concurrent radar–lidar data

3.4.1 Addressing forest disturbance or deforestation

The acquisition times between repeat-pass InSAR and spaceborne lidar usually do not overlap in practice. For example, the ALOS-2 InSAR acquisitions would match those of GEDI ideally; however, there are not many ALOS-2 InSAR pairs available in the data archive that are freely accessible. In contrast, ALOS-1 InSAR pairs are consistently available covering the study regions; however, there is an approximate 10-year time gap. To address this issue, a twofold solution is provided in Sect. 3.4. The temporal evolution of forests is primarily influenced by two factors: forest disturbance (including deforestation) and natural growth. Forest disturbance leads to forest loss and, in some cases, conversion to other land uses, such as bare ground. This change is poorly characterized by the signal model (3) and must be addressed in advance.

The height of the bare ground and short vegetation (≤10 m) (Lawrence, 2013; Allaby, 2012) is better inverted by jointly using the SAR backscatter and spaceborne lidar data compared to the InSAR correlation-based approach (Lei, 2016). Moreover, annual global backscatter products are available from the archive of JAXA, enabling the fusion of the ALOS-2 SAR backscatter information and spaceborne lidar under nearly concurrent acquisition conditions. Figure 8 presents an example of forest loss detection by SAR backscatter-based inversion within the New England region: the majority of forest disturbance areas as defined by the global forest change products (Hansen et al., 2013) were detected by the backscatter-based estimates. Based on the statistics over the New England region, 72 % of disturbed areas have been detected using the backscatter-based estimates.

Short forests undergoing rapid temporal height variations within short intervals cannot be adequately captured by current ALOS-1/2 datasets. These dynamic changes can be better resolved using dense time series data from TanDEM-X (Treuhaft et al., 2017; Lei et al., 2018), Sentinel-1 (Bhogapurapu et al., 2024), and the forthcoming NISAR mission.

Figure 11 Approximating the modified sinc model with the original model by aligning two points at 10 and 30 m: the newly fitted parameters at (a) Harvard Forest sites, (b) the forest site in Vermont, and (c) Howland Forest site. The y axis represents the coherence magnitude estimated from Eqs. (3) and (15).

Figure 12 Approximating the modified sinc model with the original model by aligning two points at either 10 and 15 m or 25 and 30 m: the newly fitted parameters at (a) the Harvard Forest site, (b) the White Mountain National Forest site, and (c) the Howland Forest site. The y axis represents the coherence magnitude estimated from Eqs. (3) and (15).

Figure 13 Illustration of the dynamic behavior of fitted (S,C) parameters for three representative forest sites when approximating the modified sinc model with original model by local fitting at fixing points at each forest height interval $\left[x,x+\mathrm{5}\right]$.

3.4.2 A modified model considering the natural forest growth

The natural growth of forests remains another concern. In temperate forest regions, the intact forest landscape (IFL), along with forest height, remains quite stable and only changes slightly as trees grow (Potapov et al., 2008; Riofrío et al., 2023). In the New England region, the growth rate was found to be inversely proportional to forest height. This finding is based on a comparison between ICESat-1 or LVIS lidar data acquired before 2009 and LVIS lidar data collected in 2021. Using these datasets, the temporal evolution of forest height at two epochs is modeled as follows:

$$\begin{array}{}\text{(11)}& h_{\mathrm{v}}\left(t_{\mathrm{2}}\right)=h_{\mathrm{v}}\left(t_{\mathrm{1}}\right)+G\cdot \left(t_{\mathrm{2}}-t_{\mathrm{1}}\right),\end{array}$$

where h_v represents the time-dependent forest height and t₁andt₂ denote the initial and subsequent epochs, respectively. Although the forest growth rate is derived by comparing pre- and post-growth states, it can be modeled based only on forest height data from either epoch, for example, using the initial forest height h_v(t₁):

$$\begin{array}{}\text{(12)}& G=a\cdot h_{\mathrm{v}}\left(t_{\mathrm{1}}\right)+b,\end{array}$$

where a and b are linear coefficients. If a dense time series of forest height data over certain forest land cover is provided, the above equation can be constructed in a differential form as

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial h_{\mathrm{v}}}{\partial t}}=a\cdot h_{\mathrm{v}}+b.\end{array}$$

This ordinary differential equation yields the following expression:

$$\begin{array}{}\text{(14)}& h_{\mathrm{v}}\left(t_{\mathrm{1}}\right)=\left({\displaystyle \frac{b}{a}}+H\left(t_{\mathrm{2}}\right)\right)e^{-a\left(t_{\mathrm{2}}-t_{\mathrm{1}}\right)}-{\displaystyle \frac{b}{a}}.\end{array}$$

This expression aligns with the Hossfeld model. However, it requires statistical data on annual growth rates, which is often unavailable on a large scale in practice. When the parameter a is small, Eq. (14) simplifies to the form presented in Eq. (12). Evaluation in the New England region indicates the absolute value of a is less than 0.02. By substituting Eq. (12) into the signal model, the following modified model is obtained:

$$\begin{array}{}\text{(15)}& \left|{\mathit{\gamma }}_{\mathrm{t}\mathrm{\&}\mathrm{v}}^{\mathrm{HV}}\left(t_{\mathrm{1}}\right)\right|=S\left(t_{\mathrm{1}}\right)\cdot \mathrm{sinc}\left({\displaystyle \frac{h_{\mathrm{v}}\left(t_{\mathrm{2}}\right)-b\cdot \mathrm{d}t}{\left(\mathrm{1}-a\cdot \mathrm{d}t\right)\cdot C\left(t_{\mathrm{1}}\right)}}\right),\end{array}$$

where $\left|{\mathit{\gamma }}_{\mathrm{t}\mathrm{\&}\mathrm{v}}^{\mathrm{HV}}\left(t_{\mathrm{1}}\right)\right|$ represents the InSAR coherence at initial time t₁, and it follows that the model is shifted and scaled with respect to original model (3); $\mathrm{d}t=t_{\mathrm{2}}-t_{\mathrm{1}}$.

As a representative example, Fig. 9 illustrates the fitted forest growth rate ($G=-\mathrm{0.0134}\cdot h_{\mathrm{v}}\left(t_{\mathrm{1}}\right)+\mathrm{0.464}$; unit: m yr⁻¹) for the forest height interval with the highest density at the Harvard Forest site in Massachusetts. This analysis compares LVIS lidar data acquired in 2009 and 2021 after filtering out disturbed forests areas ($G\le -\mathrm{0.1}$ m yr⁻¹) and short vegetation (h_v≤10 m). To illustrate how the InSAR correlation observations from ALOS-1 data are linked to GEDI forest estimates (with time gap of around 10 years), a simulation can be performed by inserting the fitted growth rate into model (15) and settings S(t₁)=0.9 and C(t₁)=11. Figure 10 presents the simulated results of original sinc model and its modified version. It follows that the forest height below b⋅dt cannot be well characterized. This value usually ranges from 8 to 15 in the investigation over the New England region. This finding also highlights the importance of utilizing backscatter information to estimate the heights of short forests in this case. We remark that natural growth functions are highly species-specific. From a practical inversion perspective, the application of the modified model requires precise detailed statistics of natural growth across various forest types on a large scale. However, such data are currently unavailable as existing spaceborne lidar datasets lack collocated measurements from two distinct time periods. In the absence of comprehensive forest growth data, this model is not yet recommended for direct large-scale use. Instead, it can be integrated into the framework of the original model by adjusting temporal parameters. The following subsection provides simulation examples to demonstrate this adaptation.

Figure 14 Density scatterplot illustration of relation between repeat-pass ALOS-1 InSAR coherence (3×10 window at 30 m grid) and GEDI RH98, and a red dashed fitted line with the slope as the pre-inversion metric.

Table 3 An example to illustrate the pre-inversion metric (slope) and their inversion performance (root mean square error, RMSE, and coefficient of determination, R²) for the available InSAR pairs in 2007 at the Howland Forest site in Maine, US.

Figure 15 30 m gridded forest height mosaic map based on ALOS-1 InSAR and GEDI RH98 metric for the New England region in the US, with a total area of 18 million ha. The color map ranges from 0 m (“blue” for bare surfaces) to 35 m (“red” for tall trees). It was projected onto the map coordinate of Shuttle Radar Topography Mission (SRTM) DEM products. The optical basemap is from © Microsoft Bing Maps.

Figure 16 30 m gridded forest height mosaic map based on the ALOS-1 InSAR and GEDI RH98 metric for the northeastern region of China, with a total area of 152 million ha. The optical basemap is from © Microsoft Bing Maps.

Figure 17 Density scatterplots comparing lidar validation data with forest height inversion estimates across multiple sites of the New England: the left panels show ALOS-1-based estimates; the right panels show ALOS-2-based estimates.

Table 4 Comparison of GLAD, ETH, and our ALOS-1-based canopy height products with airborne lidar data across all forest sites in the New England region.

Figure 18 Validation of forest height inversion results at the Howland Forest site: (a) LVIS RH98 canopy height map (30 m grid), (b) inversion extracted from ALOS-1 mosaic, and (c) ALOS-2-based inversion.

3.4.3 Approximating the modified model with the original model in the absence of natural growth data

Without detailed forest growth statistics, this subsection demonstrates that the modified sinc model can be well approximated in the framework of the original sinc model but with updated parameters $\left(S^{\prime },C^{\prime }\right)$ using simulations. This enables large-scale application achieved in the framework of the original model without detailed growth statistics.

It is almost impossible to derive S^′ and C^′ based on a direct one-to-one analytical transformation between models (15) and (3). Instead, such parameters can be determined by aligning the original model and modified model at two fixed points. For example, Fig. 11 presents the resolved parameters when the alignment is performed for two points at 10 and 30 m for global fitting at three typical forest sites (e.g., scenes with large forest height ranges). The best approximation is observed at the Howland Forest site, where the forest height changes slowly. Model mismatches are observed at 20 m at other two forest sites. This can be addressed by aligning two models for each short interval. Figure 12 presents the behavior of newly fitted parameters when aligning either a short forest interval (e.g., 10–15 m) or a tall forest interval (e.g., 25–30 m). The results suggest that the modified sinc model is better approximated through piecewise fitting, emphasizing the importance of local fitting to approximate modified models for a short height range within relatively homogeneous forest areas. Notably, S^′ increases and even exceeds 1, particularly for taller trees, deviating from its original physical definition $S^{\prime }\le \mathrm{1}$ (e.g., only due to moisture-induced dielectric decorrelation in the original model). As seen in the behavior of the fitted parameters for each short interval $\left[x,x+\mathrm{5}\right]$ in Fig. 13, the C^′ parameter is larger than the original parameter C for short vegetation but approaches C as vegetation height increases, while the S^′ parameter is close to the original S parameter for short trees and becomes biased for tall trees. For the three forest sites analyzed, the White Mountain National Forest shows a larger forest height change rate, leading to a greater deviation between the fitted parameters $\left(S^{\prime },C^{\prime }\right)$ and the original parameters (S,C).

Figure 19 The differential height maps of (a) ALOS-1-based inversion versus LVIS data and (b) ALOS-2-based inversion versus LVIS data.

3.5 The selection of proper InSAR pairs and mosaicking

Due to several decorrelation factors (induced by precipitation, human activities, etc.), the temporal decorrelation behavior over forested scenarios is complicated in the context of repeat-pass InSAR. For example, total decorrelation is possible to occur for ALOS-1 InSAR pairs with a temporal span larger than 2 months during the monsoon season. In this case, the InSAR correlation behavior is dominated by additional decorrelation sources and hence would no longer be well suited for forest height inversion.

In this context, a pre-inversion metric and a post-inversion error assessment need to be defined when multiple InSAR pairs are available to eliminate those scenes that would not work well for forest height inversion.

An example of this evaluation is illustrated in Fig. 14 and Table 3 below. Here, the fitness of temporal change model (3) for a given InSAR scene can be evaluated by testing the underlying assumption of the physical scattering model. Specifically, taller trees are more easily decorrelated over time due to the larger deviation of random motions compared to the short ones. A simple yet effective pre-evaluation metric can be attained by linear regression between the GEDI RH98 samples (only keeping forested land cover) and the coherence magnitude observations. A negative slope in the linear regression usually indicates a relevant validity of the inversion model. As shown in Fig. 14 and Table 3, InSAR pairs with a negative slope tend to yield more accurate estimates (e.g., the image pair with a slope of −4.76 gives better estimates with respect to the image pairs with slopes higher than −3). Note that while InSAR pairs with short temporal baselines may present higher correlation values, they may occasionally not be well suited for the temporal decorrelation model due to regional precipitation.

The post-inversion metric is defined using the figure of merit during the local fitting (see Eq. 10): once one pair of temporal parameters $\left(\left\{S\left(i,j\right),C\left(i,j\right)\right\}\right)$ is determined, a weighted squared summation of the differences between inverted height estimates and GEDI measurements over a regional window is given by

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& \mathit{\epsilon }\left(i,j\right)=\\ & {\displaystyle \frac{\sum _{\left(r,c\right)\mathit{ϵ}W}{\left[w(r,c)\left({\widehat{h}}_v\left(i+r,j+c\right)-h_{\mathrm{gedi}}\left(i+r,j+c\right)\right)\right]}^{\mathrm{2}}}{\sum _{\left(r,c\right)\mathit{ϵ}W}w{\left(r,c\right)}^{\mathrm{2}}}}.\end{array}\end{array}$$

It can be noted that the derived post-inversion metric ε(i,j) is also in the same coordinate system as GEDI samples. A Delaunay-triangulation-based natural neighborhood interpolation (Park et al., 2006) can be used afterwards to attain pixel-by-pixel evaluation. This assessment serves to guide the mosaicking of overlapping regions between consecutive InSAR scenes by preserving estimates with lower ε(i,j).

It should also be noted that since there are not enough InSAR pairs each year for ALOS-1 data, in this work, we only use the above-mentioned pre- and post-inversion metrics to select the best InSAR pair. However, if there are sufficient pairs from more recent spaceborne repeat-pass InSAR missions (such as 12 d for NISAR, 14 d for ALOS-4, and 4–8 d for LuTan-1), a synthetic InSAR coherence map can be generated by applying the monthly, seasonal median, or maximum operations (Kellndorfer et al., 2022).

Figure 20 (a) Global-fitting-based inversion (Lei et al., 2019) applied to the scene $\left(S=\mathrm{0.94},C=\mathrm{9}\right)$. (b) Comparing 90 m gridded maps from (a) with corresponding LVIS data. (c) Interpolated 30 m gridded GEDI height map. (d) Density scatterplot comparing the interpolated 30 m GEDI map with LVIS lidar data. (e) Density scatterplot comparing ALOS-2-based inversion results with LVIS lidar data.

Figure 21 The interpolated maps of temporal change parameters for (a) S and (b) C.

4 Results

This section begins by presenting large-scale forest height mosaic maps for the northeastern regions of the US and China and is followed by extensive validation for representative individual forest sites.

4.1 Forest height mosaic generation

The proposed inversion approach was developed as an automated open-source software serving as version 2 of the Forest Stand Height (FSH) software (https://github.com/Yanghai717/FSHv2 Lei and Yu, 2025, last access: 29 July 2025).

For generating the forest height mosaic map over the New England region (a total area of 18 million ha), over 100 ALOS-1 InSAR scenes were processed using multi-look averaging with 2 range looks and 10 azimuth looks, leading to a pixel size of 20 by 30 m, consistent with the SRTM grid. Approximately 15 million GEDI RH98 samples were used for model parameterization. The height estimates of short vegetation were replaced with the backscatter-based estimates using the ALOS-2 backscatter products in either 2019 or 2021. A few ALOS-2 InSAR pairs were also used for demonstrating radar–lidar fusion under concurrent acquisitions. Non-forest areas were masked out based on the 2021 National Land Cover Database (NLCD) products. The mosaic was projected onto the same geographic coordinate grid as the SRTM DEM product. The forest height mosaic is depicted in Fig. 15. The absence of discontinuity between adjacent scenes confirms the consistency of the forest height estimates. Additionally, the coastal region is included in the inversion despite potential challenges posed by weather conditions (as reported in Lei et al., 2019). This underscores the advantage of using the global-to-local two-stage inversion approach to handle fast spatially varying temporal change factors induced by different land covers or weather/climate conditions. Further quantitative evaluation is conducted in subsequent sections, with a focus on each individual forest site.

Figure 22 Validation of 30 m gridded forest height inversion at the Harvard Forest site: (a) LVIS lidar RH98, (b) ALOS-1-based estimates, and (c) ALOS-2-based inversion.

Figure 23 Differential height maps over the Harvard Forest site: (a) ALOS-1 mosaic versus LVIS lidar, and (b) ALOS-2 single scene versus LVIS lidar.

For the northeast of China, 688 ALOS-1 InSAR scenes and 160 million GEDI samples were used to generate the mosaic covering the 5 provinces (total area of 152 million ha). Non-forested areas were masked out based on the 2021 ESRI global land cover maps. ALOS-2 global backscatter maps for 2019 and 2020 were employed to estimate the height of short trees to match the acquisition time of the validation airborne lidar data. The final forest height mosaic is shown in Fig. 16. It is noted the area outside the coverage of GEDI observation (> 51.6° N) was discarded.

The generated products are made available via https://doi.org/10.5281/zenodo.11640299 (Yu and Lei, 2024). Further evaluation is shown in Sect. 4.2 and 4.3. In both cases, small values of κ_z are maintained for all available InSAR pairs (κ_z are below 0.15 rad m⁻¹, and the mean values are 0.032 and 0.029 rad m⁻¹ for Chinese and American datasets, respectively), which conforms to the assumption made in Lei and Siqueira (2014).

4.2 Validation over the New England region

This subsection presents the validation of forest height inversion across representative test sites in New England, US. First, we assess the accuracy of the inversion results for all selected forest sites using density scatterplots and corresponding error metrics including the root mean square error (RMSE), coefficient of determination (R²), standard deviation (SD), and bias. Unless otherwise stated, all density scatterplots and associated error metrics in this study are based on a 0.8 ha (3×3 pixel) aggregated pixel size. A comparative analysis is then conducted between our inversion results and two existing forest height products: (1) the 30 m resolution GLAD canopy height map, derived by fusing GEDI and Landsat time series (Potapov et al., 2021), and (2) the 10 m resolution ETH canopy height map, generated by fusing GEDI and Sentinel-2 data (Lang et al., 2022). Both products are validated against reference lidar data.

Subsequently, case studies are presented to evaluate inversion performance under varying conditions. The Howland Forest site is analyzed to quantify improvements over our earlier methodology. A second case study focuses on the high-biomass region at the Harvard Forest site to assess inversion robustness in dense canopies. Finally, validation is extended to the White Mountain National Forest (WMNF) site with hilly topography using high-resolution small-footprint lidar data as the validation data.

Figure 24 (a) An example of a histogram formed by small-footprint CHM values within the GEDI footprint, with the mean height and 98th percentile height marked by the red dashed and red solid lines, respectively. (b) 30 m gridded reprocessed forest height based on the ERH50 metric, (c) 30 m gridded reprocessed forest height based on the ERH98 metric, and (d) corresponding forest height estimates extracted from ALOS-1 mosaic.

Figure 25 Differential height map between the ALOS-1 mosaic and GRANIT ERH98 map.

4.2.1 Summary of the validation results over the northeastern US

Figure 17 presents density scatterplots comparing forest height estimates (derived from ALOS-1 mosaics or ALOS-2 single-scene inversions where suitable InSAR pairs are available) across all test sites in the New England region of the US. Error metrics for these estimates are reported in the corresponding scatterplots.

In short, the proposed inversion approach is capable of estimating forest height with an RMSE of 3–4 m in areas such as Howland and Harvard Forest sites, characterized by relatively flat topography and minimal human activity influence. In contrast, hilly or suburban areas, such as WMNF, GMNF, and Naugatuck State Forest (NSF) sites, exhibit slightly lower accuracy (RMSE 4–5 m). The ALOS-2-based inversion generally presents superior performance due to enhanced InSAR correlation behavior resulting from shorter temporal baselines and less temporal discrepancy between radar and lidar data. For ALOS-1 inversion, the time mismatch between ALOS-1 and GEDI data is not a fatal problem as the inversion is carried out for temperate regions where the intact forest landscapes and forest height (of mature forests) remain stable. While our solution for addressing forest growth may not fully resolve temporal uncertainties, the resulting errors remain relatively minor, as evidenced by an RMSE of 3–4 m. This is also supported by the finding that the ALOS-1-based inversion is occasionally more accurate than ALOS-2-based estimates.

In addition, we evaluate the inversion performance of two widely recognized global forest products: the GEDI-Sentinel (ETH) product (Lang et al., 2022) and the GEDI-Landsat (GLAD) product (Potapov et al., 2021). As shown in Table 4, both products exhibit significant biases across forest sites in the New England region compared with our ALOS-1-based estimates. Specifically, the GLAD product systematically underestimates canopy height, likely due to saturation effects in dense canopies, while the ETH product demonstrates larger systematic biases, consistently overestimating canopy height at all sites. These findings are consistent with the analysis by Qi et al. (2025). In contrast, our inversion method achieves lower RMSE and smaller biases in most cases, with the exception of Naugatuck State Forest. At this site, hilly topography and suburban land cover likely contribute to reduced inversion accuracy. Notably, the ETH product exhibits lower SD and higher height-related R² values, which is likely attributable to its integration of high-resolution Sentinel-2 data.

Figure 26 Density scatterplots comparing lidar validation data with forest height inversion estimates across multiple sites of northeastern China based on the ALOS-1 InSAR observation (left panels) and the ALOS-2 observation (right panels).

4.2.2 Howland Forest site

The Howland Forest was selected as one of the representative test sites, continuing from previous efforts in developing the inversion approaches (Lei et al., 2018; Lei and Siqueira, 2022). Comparing results from these earlier studies with the current work allows for an evaluation of performance improvements. A strip of LVIS lidar data acquired in 2009 was used as the reference to assess the inversion accuracy of both the ALOS-1-based forest height mosaic map and the ALOS-2-based inversion from a single pair of ALOS-2 data (frame: 890, orbit: 37). The comparison between the inverted forest height estimates, and the LVIS lidar data is presented in Fig. 18, generating the differential height maps between the inversion and the lidar data as shown in Fig. 19a and b. For quantitative analysis, density scatterplots and corresponding statistical error metrics comparing inversion results with validation data are displayed in the first row of Fig. 17.

Figure 27 Comparison of the forest height inversion with lidar data at the Hubao National Park site: (a) ALS ERH98 metric map. The red rectangle denotes the coverage of (b) ALOS-1-based single-scene inversion, whereas the blue rectangle indicates the coverage of (c) ALOS-2-based single-scene inversion.

Figure 28 Differential height map between (a) the ALOS-1-based inversion and ERH98 validation data and between (b) the ALOS-2-based inversion and ERH98 validation data.

The comparison reveals that the ALOS-1-based mosaic inversion estimates forest height with an RMSE of 3.8 m. In contrast, the ALOS-2-based single-scene inversion achieves enhanced accuracy (RMSE = 3.6 m), likely due to improved correlation from its 14 d temporal baseline. To compare these results against our earlier work (Lei et al., 2019), we applied the global-fitting-based inversion method (in the left panel of Fig. 20), which yielded an inversion accuracy of RMSE = 4.38 m at an aggregated pixel size of 0.81 ha. This demonstrates that the global-to-local two-stage inversion approach significantly improves inversion accuracy, particularly over tall forest regions.

By synergistically combining SAR data with GEDI samples, this methodology not only resolves the wall-to-wall mapping limitations inherent to GEDI's discrete sampling, but also improves inversion precision. As shown in the right panel of Fig. 20, the 30 m interpolated GEDI-based forest height maps still face discontinuity problems. However, an accuracy improvement of up to 20 % has been achieved for ALOS-2-based estimates. This enhancement stems from the refined characterization of temporal parameters (S,C) enabled by leveraging GEDI's dense spatial sampling, as shown in Fig. 21. As anticipated in Sect. 3.4.3, the saturation behavior of S parameters arises from approximating modified forest growth model (15) using original model (3) in the local window with taller forests.

Figure 29 Comparison of the forest height mapping over the Genhe Forest Bureau: (a) ALS ERH98 metric generated based on 25 m footprint, (b) ALOS-1 mosaic with red rectangle box indicating the overlapping area with (a), and (c) ALOS-2-based single-scene inversion.

Figure 30 Differential height map between (a) the ALOS-1-based inversion and ERH98 validation data and between (b) the ALOS-2-based inversion and ERH98 validation data.

4.2.3 Harvard Forest site

The Harvard Forest site was selected to evaluate the inversion in a region characterized by high biomass up to 400 Mg ha⁻¹ (Tang et al., 2021). Figure 22 presents the LVIS validation data acquired in 2021 covering the Harvard Forest area in panel (a) and the forest height estimates extracted from the ALOS-1 mosaic and the ALOS-2 single-scene inversion results (frame: 2770, orbit: 141) in panels (b) and (c). The comparison of ALOS-1 and ALOS-2 inversion results against the validation data is illustrated in the differential height maps in Figure 23. The density scatterplots from these two comparisons are given in the second row of Fig. 17.

Both the ALOS-1 mosaic and ALOS-2 single-scene inversion are capable of estimating forest height with an RMSE of 4 m. The biased estimates occurred in taller forest stands may be attributed to the degraded sensitivity of GEDI measurements over dense forest stands (Fayad et al., 2022).

4.2.4 White Mountain National Forest site

The evaluation of the forest height inversion will also be extended to mountainous areas such as the WMNF site, considering the potential challenges GEDI and InSAR observations might face in these regions. These challenges include GEDI's geolocation shifts and slope effects as well as the radar's viewing geometry problems (e.g., layover, shadow, and foreshortening).

A high-resolution canopy height model (CHM), derived from small-footprint GRANIT lidar data, serves as the validation reference. As outlined in Sect. 2.2.3, due to the footprint difference between small-footprint lidar data and GEDI observations, an equivalent RH metric must be extracted within the same footprint size as GEDI observations to ensure comparability. Without this adjustment, directly comparing GEDI-based forest height inversions with the reprojected CHM (via resampling or multi-pixel averaging) introduces significant bias. Figure 24 illustrates the difference between the ERH98 metric and the mean value (or ERH50 metric) in panel (a) and shows the maps of these two metrics in panels (b) and (c). Notably, ALOS-1-based forest height estimates (Fig. 24d) align closely with the ERH98 metric (Fig. 24b). The differential map between the ALOS-1 mosaic and GRANIT ERH98 map is shown in Fig. 25, and the associated density scatterplot can be found in the third row of Fig. 17.

4.3 Validation against ALS data over the northeastern region of China

The forest height mosaic for northeastern China was validated exclusively against small-footprint airborne laser scanning (ALS) data at representative forest sites. These high-resolution ALS datasets were processed into ERH98 metric maps following the methodology detailed in Sect. 2.2.3. Initial analysis focused on density scatterplots across all surveyed sites. For a deeper investigation, two case studies were examined: Hubao National Park and Genhe Forest Bureau. Hubao National Park was selected due to the absence of significant forest disturbance, as demonstrated in Fig. 4, while Genhe Forest Bureau was chosen for testing the inversion performance over the boreal forest bioregion in China.

4.3.1 Summary of the validation over the northeast of China

Forest height estimates across all sites across the northeast of China were validated against ERH98 metrics derived from small-footprint airborne lidar data. Density scatterplots as well as error metrics for the representative forest sites are summarized in Fig. 26. The highest accuracy was observed at the Mengjiagang Forest site, achieving an RMSE of 3.32 m and an R² of up to 0.84. Most inversions exhibit a slight negative bias, likely due to GEDI's reduced signal penetration capability compared to airborne lidar. Slightly less accurate estimates are provided by the ALOS-2-based inversion at the Saihanba forest site, which is attributed to the limited overlapping area between the ALOS-2 single-scene inversion and ALS validation observations. This limitation arises from the distribution of heterogeneous land cover influenced by human activities. Overall, forest height estimates align closely with ERH98 lidar benchmarks, with accuracies of 3–4 m (even below 3.5 m at three sites) and R² predominantly exceeding 0.65.

Notably, inversion performance in northeastern China surpasses results from GEDI calibration sites in the northeastern US, likely due to fewer forest disturbance in the selected represented forest sites (See Fig. 4) compared to those in the northeastern US, as shown in Fig. 3.

4.3.2 Hubao National Park Forest site

Hubao National Park is selected as it represents one of the typical temperate regions with the richest biodiversity in terms of wildlife and plants in the Northern Hemisphere. Figure 27a displays an ERH98 map derived from reprocessed 1 m resolution canopy height model (CHM) data acquired in 2018 using a window size same as GEDI footprints. Panels (b) and (c) present forest height estimates from ALOS-1 single-scene inversion (Frame 860, Orbit 421) and ALOS-2 single-pair InSAR data (Frame 860, Orbit 130), respectively. The detailed differential height maps are provided in Fig. 28.

As shown in the density scatterplots (see the bottom row of Fig. 26), both ALOS-1 mosaic and ALOS-2 single-scene inversions align closely with ALS-derived ERH98 data, achieving accuracies of 3.5–3.8 m with an R² of up to 0.7. These results demonstrate the inversion precision comparable to airborne lidar measurements across both short and tall vegetation types. Differential height maps highlight discrepancies in transitional zones between forested and bare surfaces, underscoring the need for higher-spatial-resolution data integration (e.g., fusing TDX-GEDI data; Hu et al., 2024; Lei et al., 2021; Qi et al., 2025) to refine estimates in such areas.

Slightly better performance for ALOS-1-based inversion is attributed to the fact that the ALOS-1 InSAR data archive offers the possibility of picking out the best InSAR pair with better correlation behavior; however, the availability of proper ALOS-2 InSAR data is limited.

4.3.3 Genhe Forest Bureau

A second case study examines the Genhe Forest Bureau situated in northeastern China within one of the country's northernmost boreal forest regions. Figure 29 presents the (a) reference ERH98 map derived from 1 m resolution airborne lidar data, (b) forest height estimates from the ALOS-1-based mosaic, and (c) results from the ALOS-2 single-scene inversion. Differential height maps are shown in Fig. 30. The spatial distribution of forest heights in both ALOS-1 and ALOS-2 inversions closely matches the reference ERH98 map, though primary discrepancies occur in short vegetation and bare surfaces, which is likely linked to anthropogenic disturbances.

Quantitative validation (see the third row of Fig. 26) confirms strong agreement between inverted and reference forest height estimates, with an RMSE of 3.6 m and an R² of 0.65, demonstrating the method's effectiveness and accuracy in boreal ecosystems.

5 Discussion

This section outlines key limitations of the proposed inversion in the framework. A primary challenge stems from the complex InSAR correlation behavior induced by weather fluctuations (e.g., precipitation), particularly when constrained by a limited number of InSAR pairs. Future missions, such as NISAR and BIOMASS (Quegan et al., 2019), are expected to mitigate this issue. As data stacks accumulate rapidly within each season, seasonally synthesized coherence maps could be produced by averaging or selecting maximum coherence values from all available pairs (Kellndorfer et al., 2022). For example, leveraging NISAR's 12 d repeat cycle would generate 7 interferometric pairs per season (or 30 annually), enabling the creation of seasonally or annually averaged coherence maps. These refined datasets could significantly enhance the robustness of forest height inversion.

To address the temporal mismatch between ALOS-1 and GEDI datasets, we propose a twofold solution to account for forest change dynamics. While this approach does not fully capture inherent forest variability, the achieved forest height estimation accuracy of 3–4 m ha⁻¹ suggests that its impact is negligible in temperate regions, where intact forest landscapes exhibit minimal canopy height variation due to gradual tree growth.

The methodology also depends on precise, spatially representative calibration samples from GEDI. Two key challenges arise: (1) slope-induced biases in GEDI forest height estimates and (2) sparse sampling coverage in boreal and equatorial tropical regions. The first limitation can be mitigated using the RH metric derived from slope-corrected waveforms (Wang et al., 2019). The second issue may be resolved by integrating complementary lidar datasets, such as NASA's ICESat-1/2 missions. Combining ALOS-1/2 data with both GEDI and ICESat-1/2 observations would improve the method's accuracy and adaptability, particularly in tropical ecosystems.

Finally, the current inversion framework employs a fixed-size, distance-based weighting window to perform local fitting. This approach could be enhanced using an adaptively spatially varying window size driven by multi-parameter classification.

6 Data availability

The forest height mosaics over the northeastern parts of US and China are available at https://doi.org/10.5281/zenodo.11640299 (Yu and Lei, 2024). The used ALOS-1 data can be found via the Alaska Satellite Facility at https://search.asf.alaska.edu/#/?dataset=_ALOS (Alaska Satellite Facility, 2025). GEDI data (from 2018 to 2023) can be downloaded from the EARTHDATA SEARCH website at https://search.earthdata.nasa.gov/ (National Aeronautics and Space Administration, 2025a). Regarding the validation data, the LVIS and GRANIT lidar data can be found at https://lvis.gsfc.nasa.gov/Home/index.html (National Aeronautics and Space Administration, 2025b) and https://lidar.unh.edu/map/, University of Hampshire, 2025.

7 Code availability

The forest height mosaics are generated using the following software tools. First, ISCE Version 2.4+ (https://github.com/isce-framework/isce2, Jet Propulsion Laboratory, 2025; in particular the stripmapApp function) is used to preprocess the two ALOS-1/-2 images for producing geocoded interferometric coherence maps. Then, FSH software Version 2 (https://doi.org/10.5281/zenodo.16547783, Lei and Yu, 2025) is used to invert forest height by fusing GEDI and InSAR data and perform the mosaicking. Several preprocessing steps utilize basic Python libraries from FSH Version 1: https://github.com/leiyangleon/FSH; leiyangleon, 2025.

8 Conclusions

This paper presents a global-to-local two-stage forest height inversion approach for large-scale forest stand height mapping using L-band spaceborne repeat-pass InSAR and spaceborne GEDI lidar. This work extended our previous efforts in forest stand height mapping (FSH: https://github.com/leiyangleon/FSH, last access: 29 July 2025; Lei and Siqueira, 2014; Lei and Siqueira, 2015; Lei et al., 2018) at a large scale by incorporating GEDI lidar samples for capturing local information. The sparsely yet extensively distributed lidar samples provided by the GEDI mission are used to parameterize the semi-empirical InSAR scattering model and to obtain forest height estimates. Building on earlier works (Lei et al., 2019; Lei and Siqueira, 2022), this paper removed the assumptions previously imposed by the limited availability of calibration samples before and developed a new inversion approach based on a global-to-local two-stage inversion scheme. An effective use of regional GEDI samples in this approach allows for a finer characterization of temporal decorrelation patterns and thus higher accuracy in forest height inversion while also suppressing issues in individual GEDI sample, e.g., geolocation errors. This approach was supported by fusing the ALOS-2 InSAR data and GEDI data under nearly concurrent conditions. This approach is further applied to open-access ALOS-1 data for testing its mapping capabilities at a large scale. To address the temporal mismatch between ALOS-1 data and GEDI data, the introduction of fusing ALOS-2 backscatter data and GEDI data provides us with the ability to detect disturbed forest areas. Furthermore, a modified signal model is derived for addressing natural forest growth over temperate forest regions, where the intact forest landscape and forest height are stable with slight changes. In the absence of detailed forest growth data, simulations confirmed that the modified model can be effectively approximated by the original model's form through local fitting. For evaluating its performance, two forest height mosaic maps were generated to investigate the northeastern regions in the US and China, covering a total area of 18 and 152 million ha, respectively. Validation of the forest height estimates demonstrates substantial accuracy improvements achieved by the proposed approach compared to the previous efforts, i.e., from 4 m RMSE for 6.25 ha aggregated pixel size to 3.8 m RMSE for 0.81 ha aggregated pixel size at the Howland Forest site. The proposed fusion approach not only addresses the sparse spatial sampling problem of the GEDI mission, but also improves the accuracy of forest height estimates compared to the GEDI interpolated height estimates by 20 %. The extensive evaluation of forest height inversion against LVIS lidar data over the northeastern US indicates an accuracy of 3–4 m on the order of 0.81 ha over smooth areas and 4–5 m over hilly areas, while the forest height estimates over northeastern China compare well with small-footprint lidar validation data even at an accuracy of below 3.5 m on the order of sub-hectares and with an R² mostly above 0.6.

Despite the limitations outlined in Sect. 5, the method demonstrates promising forest height accuracy at the sub-hectare scale using open-access InSAR and lidar datasets, underscoring its potential for integration with current and future missions (e.g., ALOS-4, NISAR, LuTan-1, MOLI, TECIS). This framework offers a cost-effective solution for large-scale forest height mapping based on open-access remote sensing data, particularly in regions lacking multi-baseline bistatic InSAR data for advanced techniques such as PolInSAR or TomoSAR.