Improved Model-Based Forest Height Inversion Using Airborne L-Band Repeat-Pass Dual-Baseline Pol-InSAR Data

: This paper proposes an improved model-based forest height inversion method for airborne L-band dual-baseline repeat-pass polarimetric synthetic aperture radar interferometry (PolInSAR) collections. A two-layer physical model with various volumetric scattering attenuation and dynamic motion properties is ﬁrst designed based on the traditional Random Motion over Ground (RMoG) model. Related PolInSAR coherence functions with both volumetric and temporal decorrelations incorporated are derived, where the impacts of homogenous and heterogeneous attenuation and dynamic motion properties on the performance of forest height inversion were investigated by the Linear Volume Attenuation (LVA), Quadratic Volume Attenuation (QVA), Linear Volume Motion (LVM), and Quadratic Volume Motion (QVM) depictions in the volume layer. Dual-baseline PolInSAR data were acquired to increase the degree of freedom (DOF) of the coherence observations and thereby provide extra constraints on the forest parameters to address the underdetermined problem. The experiments were carried out on a boreal forest in Canada and a tropical one in Gabon, where physical models with LVA + QVM (RMSE: 3.56 m) and QVA + LVM (RMSE: 6.83 m) exhibited better performances on the forest height inversion over the boreal and tropical forest sites, respectively. To leverage the advantages of LVA, QVA, LVM, and QVM, a pixel-wise optimization strategy was used to obtain the best forest height inversion performance for the range of attenuation and motion proﬁles considered. This pixel-wise optimization surpasses the best-performing single model and achieves forest height inversion results with an RMSE of 3.21 m in the boreal forest site and an RMSE of 6.48 m in the tropical forest site.


Introduction
Forest height has been identified as a key biophysical parameter to quantify the forest Above-Ground Biomass (AGB), which can be further used to justify the role of forests in the global carbon cycle and ecosystem productivity. Remote sensing techniques like Light Detection And Ranging (LiDAR) and Polarimetric Interferometric Synthetic Aperture Radar (PolInSAR) have been widely applied for forest height estimation in the past decades. LiDAR measurements are normally with higher vertical accuracies, but the availability is limited by the acquisition cost and persistent cloud cover. Compared to LiDAR, PolInSAR measurements are normally with better spatial consistency over a large scale, but the vertical accuracies largely rely on the model used in the inversion process. Even so, we still consider PolInSAR as a powerful tool for forest height estimation as related measurements are sensitive to both physical properties and vertical structures of the forest medium [1]. Prior studies estimated the canopy height through a series of model-based methods applied to PolInSAR data. These sophisticated physical models establish mathematical relations between the coherence observations and forest parameters from which the forest height can Remote Sens. 2022, 14, 5234 2 of 20 be estimated through the solution of non-linear equations [2][3][4][5][6]. One most used model, the Random Volume over Ground (RVoG) model [7], represents the forest as a homogeneous volume of randomly oriented scattering particles statically distributed over the ground surface. In this model, the volumetric attenuation of radar scattering amplitude follows an exponential distribution with a constant extinction coefficient, and the derived PolInSAR coherence function is fully determined by four parameters, including ground phase, forest height, mean extinction, and ground-to-volume amplitude ratio. Several complementary models have also been devised based on the RVoG model by integrating the vertical heterogeneity into the forest representation [3,8]. These models have been applied to various forest types and sensor wavelengths using single-pass PolInSAR data [5,[9][10][11][12][13].
In the repeat-pass interferometry, the temporal decorrelation derived from the dynamic changes in the vegetation and ground properties occurring between acquisition times affects the performance of forest parameter inversion. Therefore, subsequent studies focused on overcoming the impacts of temporal decorrelation through quantifying and compensating temporal components in the repeat-pass PolInSAR coherence. The volumetric temporal decorrelation (RVoG + VTD) model was first proposed in [14], describing the temporal decorrelation of volume and ground layers as two complex coefficients embedded in the coherence function of the RVoG model. More physically based models, such as the Random Motion over Ground (RMoG) model [15] were designed to compensate for the temporal decorrelation components derived from the scatterer motions. A physical model accounting for both the position and dielectric property changes [2] was proposed based on the RMoG model, where the dielectric property changes of the volume and ground layers were depicted by two complex coefficients integrated into the RMoG coherence function.
Each physical model mentioned above has its unique advantages when dealing with different types of PolInSAR data. The data involved in this study were collected using Uninhabited Aerial Vehicle Synthetic Aperture Radar (UAVSAR) [16], which is a fully polarimetric L-band radar system collecting repeat-pass interferometric data within short temporal baselines (half an hour to hours) under various weather conditions. When the data were collected with suitable weather conditions, the dielectric changes derived from the moisture content variations in both the canopy and ground layers were negligible, and it was reasonable to assume a stationary ground in the model. In this condition, the primary temporal decorrelation was assumed to be derived from the wind-induced movements of the scatterers in volume layers, as in the RMoG model [17]. Though this model has been applied to the forest height estimation using UAVSAR data in several studies [15,18], the effect of different scattering attenuation and random motion characteristics on the inversion accuracy has not been explored yet. This is especially relevant when dealing with Pol-InSAR data acquired at lower frequencies, such as the L-band and P-band, where longer radar waves are more likely to interact with the scattering elements within the canopy and, thereby, be more sensitive to the forest structure [3]. To this end, we considered both homogenous and heterogeneous scattering attenuation and random motion properties (LAV, QVA, LVM, QVM) in the proposed model and investigated their performance on the forest height inversion over different forest sites. In addition to that, a pixel-wise optimization strategy was further used to leverage the advantages of different attenuation and motion profiles.
The large footprint full-waveform Light Detection and Ranging (Lidar) data [19,20] obtained by airborne Land, Vegetation, and Ice Sensor (LVIS) over the two forest sites were used as the ground truth for the result analysis. This has been justified in several studies through the cross-validation between LVIS height metrics and the field data [21,22]. Though these LiDAR observations are sparser than the spatial resolution and coverage of radar images, they can characterize a variety of vegetation and terrain types with high vertical accuracy. Therefore, down-sampling processes are conducted on the retrieved forest height for a better comparison during the result analysis, which will be indicated in detail in the following. This paper is organized as follows. The test sites and datasets used in the experiments are shown in Section 2. Section 3 describes the proposed physical model, PolInSAR coherences, and forest height inversion method. Forest height inversion results over the boreal and tropical forest sites are presented in Section 4, followed by a discussion in Section 5. Finally, conclusions are drawn in Section 6.

Test Sites
The test sites covered by boreal and tropical forests were involved in the experiment to investigate the applicability of the physical model with different volumetric profiles. The boreal forest site is part of the Boreal Ecosystem Research and Monitoring Sites (BERMS) in the Saskatchewan Province of Canada, covering around 16 km along the east-west direction (104 • 46 38 W to 105 • 1 15 W) and 9 km along the north-south (53 • 52 29 N to 53 • 56 10 N). The southwest corner of the forest site (Candle Lake) is located about 70 km northeast of Prince Albert. Dominant tree species at the BERMS site are jack pine, black spruce, and aspen, with mean canopy heights of around 14 m, 11 m, and 21 m [23,24]. The topography is relatively smooth, with a local elevation ranging from 504 m to 592 m above mean sea level (AMSL), and the mean annual precipitation is about 466 mm [25,26].
The tropical test site is in the Pongara National Park (PNP) on the south bank of the Gabon Estuary. The image dataset covers a total area of about 12,008 ha within 0 • 4 15 N to 0 • 7 50 N and 9 • 45 20 E to 9 • 55 1 E. The dominant tree species in the PNP site are red mangroves covering about 55% of the total area [27]. Growing in brackish water with a lot of sediment, some of the mangroves tower up to 60 m. The PNP site also exhibits a relatively smooth topography, with an elevation range of less than 104 m. The annual precipitation over this area ranges from 2400 mm on the eastern portion to 2830 mm on the western side [28], and the annual mean temperature ranges in a small interval averaging approximately 26.3 ± 0.9 • C [29].

Datasets
Three types of data (airborne SAR, airborne full-waveform LiDAR, and simulated SAR data) were involved in developing the forest height inversion algorithm. The simulated SAR data were employed to evaluate the sensitivity of the different attenuation profiles in the proposed model, and related airborne acquisitions were used as the real data for the validation of the proposed forest height inversion algorithm. The airborne data over the boreal forest site were collected within the Arctic Boreal Vulnerability Experiment (ABoVE) campaign deployed by NASA [30]. As part of the ABoVE collection, the airborne SAR data over the BERMS site were captured in late August 2018, and the airborne LiDAR data were collected by LVIS in 2019. The airborne data over the tropical site were acquired within the AfriSAR campaign in 2016 [19], which was deployed in support of NASA-ISRO Synthetic Aperture Radar (NISAR), Global Ecosystem Dynamics Initiative (GEDI), and the ESA's BIOMASS mission [31][32][33]. The related SAR and LiDAR data used in this study were both collected over PNP in early February 2016.

UAVSAR
UAVSAR is a fully polarimetric L-band (1.26 GHz, 80 MHz bandwidth) SAR system deployed and operated by NASA's Jet Propulsion Laboratory (JPL) [16]. The spatial resolution of the UAVSAR single-look complex (SLC) image is 0.6 m in the azimuth direction and 1.6 m in the slant range, which is averaged by a rectangular window with sizes of eight pixels in azimuth and two pixels in the slant range to reduce speckle noise and to generate multi-looked products with pixel spacings of 4.8 m in azimuth and 3.2 m in slant range. The UAVSAR data over the BERMS site were collected in eight tracks with uniformly distributed vertical separations of

LVIS
LVIS is a medium-altitude imaging laser altimeter that provides accurate vertical measurements of the canopy top and underlying topography by digitally recording the returned signals [20]. Relative Height 100 (RH100) metrics in the LVIS Level-2 collection are typically associated with maximum tree height within a resolution beam [21,22]. This study used RH100 metrics as the ground truth to evaluate the performance of forest height estimation. RH100 has been demonstrated to be a suitable forest height reference with RMSEs of 3.29 m in boreal forests and 2.87 m in tropical forests when compared with field data [21,22] spacings providing temporal baselines of 25-175 min between the pairs. Five repeat-pa tracks were acquired over the PNP site with nonuniform vertical baselines (0, 20, 45, a 105 m) and roughly uniform temporal baselines from 12:30:11 to 14:13:21 UTC in 25 m spacings on 27 February 2016. The geocoded multi-looked SAR images in the Pauli co combination over BERMS and PNP are presented in Figures 1a and 2a, respectively.

LVIS
LVIS is a medium-altitude imaging laser altimeter that provides accurate verti measurements of the canopy top and underlying topography by digitally recording t returned signals [20]. Relative Height 100 (RH100) metrics in the LVIS Level-2 collecti are typically associated with maximum tree height within a resolution beam [21,22]. Th study used RH100 metrics as the ground truth to evaluate the performance of fore height estimation. RH100 has been demonstrated to be a suitable forest height referen with RMSEs of 3.29 m in boreal forests and 2.87 m in tropical forests when compar with field data [21,22]. The LVIS data over BERMS (10 m footprint diameter) were c lected in July-August 2019 with a flight altitude of 12.5 km and a swath width of 2.

PolSARproSim+
The simulated SAR data generated by PolSARproSim+ were used to explore the scattering attenuation properties of the different tree species. PolSARproSim+ can be considered an improved version of the ESA's PolSARproSim software, which allows more flexible simulations of different mission requirements, observational scenes, and instrument parameters [34,35]. For example, this package allows the simulation of mixed forest stands with an external definition of allometric parameters while introducing temporal decorrelation effects originating from the rain and wind. These extensions provide a powerful tool for investigating the scattering attenuation properties inside the forest and temporal decorrelation effects on polarimetric interferometry and tomography. In this study, SAR simulations were performed using the system parameter of UAVSAR over 14 deciduous and four coniferous forest stands, as listed in Table 1.

PolSARproSim+
The simulated SAR data generated by PolSARproSim+ were used to explore the scattering attenuation properties of the different tree species. PolSARproSim+ can be considered an improved version of the ESA's PolSARproSim software, which allows more flexible simulations of different mission requirements, observational scenes, and instrument parameters [34,35]. For example, this package allows the simulation of mixed forest stands with an external definition of allometric parameters while introducing temporal decorrelation effects originating from the rain and wind. These extensions provide a powerful tool for investigating the scattering attenuation properties inside the forest and temporal decorrelation effects on polarimetric interferometry and tomography. In this study, SAR simulations were performed using the system parameter of UAVSAR over 14 deciduous and four coniferous forest stands, as listed in Table 1.

Physical Model
As shown in Figure 3, the forest stands were simulated as volumes of dynamic scatterers with either homogeneous or heterogeneous attenuation and random motion over a stable ground in the condition of repeat-pass UAVSAR PolInSAR configuration. The attenuation is due to the propagation of waves into and out of the scattering medium in the distorted Born approximation [36], and the random motion is caused by wind-derived location and orientation perturbations [18]. The dynamic volume with homogeneous scattering attenuation and random motion has been previously modeled in the RMoG model [15]. Dynamic volumes with heterogeneous attenuation and motion are given in Figure 3, depicted by the gradient color in the volume and arrows with different lengths.

Physical Model
As shown in Figure 3, the forest stands were simulated as volumes of dynamic scatterers with either homogeneous or heterogeneous attenuation and random motion over a stable ground in the condition of repeat-pass UAVSAR PolInSAR configuration. The attenuation is due to the propagation of waves into and out of the scattering medium in the distorted Born approximation [36], and the random motion is caused by windderived location and orientation perturbations [18]. The dynamic volume with homogeneous scattering attenuation and random motion has been previously modeled in the RMoG model [15]. Dynamic volumes with heterogeneous attenuation and motion are given in Figure 3, depicted by the gradient color in the volume and arrows with different lengths. In the RMoG model, homogeneous scattering attenuation (LVA) is described by an exponential function with a constant mean extinction, and the homogeneous random motion of the dynamic scatterers is integrated as a Gaussian-statistic function with zero mean and vertically linear variance (LVM). The constant mean extinction and the vertically linear variance (constant gradient) are assumed based on characteristics of the homogeneous volume, where the scattering particles are uniformly distributed and ran- In the RMoG model, homogeneous scattering attenuation (LVA) is described by an exponential function with a constant mean extinction, and the homogeneous random motion of the dynamic scatterers is integrated as a Gaussian-statistic function with zero mean and vertically linear variance (LVM). The constant mean extinction and the vertically linear variance (constant gradient) are assumed based on characteristics of the homogeneous volume, where the scattering particles are uniformly distributed and randomly oriented. In this condition, the scattering attenuation ρ v and random motion η v profiles of the volume layer are given as follows: where σ LVA is the constant mean extinction, θ is the incidence angle, h v represents the forest height, and τ LV M is the gradient of the Gaussian motion variance. The suitability of a different functional form for describing the attenuation and motion profiles is also explored in this paper.
where σ QVA is the gradient of vertically linear extinction and τ QV M is the second derivative of the Gaussian motion variance. Equations (3) and (4) provides a heterogeneous scattering attenuation in the form of an exponential function with a vertically linear extinction (QVA) and a heterogeneous random motion in the form of Gaussian distribution with zero mean and vertically quadratic variance (QVM). The vertically linear extinction and the vertically quadratic variance (linear gradient) are provided based on the gradient structure of the volume layer in height, as depicted in Figure 3.
The scattering attenuation and random motion shape evolutions in (1)-(4) are plotted versus volume heights ranging from 0 m to 25 m for an incidence angle of 40 • in Figure 4.
domly oriented. In this condition, the scattering attenuation and random motion profiles of the volume layer are given as follows: where is the constant mean extinction, is the incidence angle, ℎ represents the forest height, and is the gradient of the Gaussian motion variance. The suitability of a different functional form for describing the attenuation and motion profiles is also explored in this paper.
where is the gradient of vertically linear extinction and is the second derivative of the Gaussian motion variance. Equations (3) and (4) provides a heterogeneous scattering attenuation in the form of an exponential function with a vertically linear extinction (QVA) and a heterogeneous random motion in the form of Gaussian distribution with zero mean and vertically quadratic variance (QVM). The vertically linear extinction and the vertically quadratic variance (linear gradient) are provided based on the gradient structure of the volume layer in height, as depicted in Figure 3.
The scattering attenuation and random motion shape evolutions in (1)-(4) are plotted versus volume heights ranging from 0 m to 25 m for an incidence angle of 40° in Figure 4.

PolInSAR Coherence Function
As an evaluation of the interferometric quality, the complex PolInSAR coherence is statistically calculated as follows: where ⃗ ⃗ 1 and ⃗ ⃗ 2 are complex unitary vectors in a certain transmit and receive polarization, 12 is the 3 × 3 non-Hermitian complex matrix, including both interferometric and polarimetric information, 11 and 22 denotes the 3 × 3 Hermitian coherency matrices describing the polarimetric properties, ̂ is the expected interferometric phase, and |̂| rep-

PolInSAR Coherence Function
As an evaluation of the interferometric quality, the complex PolInSAR coherence is statistically calculated as follows: where → ω 1 and → ω 2 are complex unitary vectors in a certain transmit and receive polarization, Ω 12 is the 3 × 3 non-Hermitian complex matrix, including both interferometric and polarimetric information, T 11 and T 22 denotes the 3 × 3 Hermitian coherency matrices describing the polarimetric properties,φ is the expected interferometric phase, and |γ| represents the coherence magnitude related to the phase noise [37]. After compensating the spectral and signal-to-noise ratio (SNR) decorrelations [12], the volumetric and temporal decorrelations become two main components in repeat-pass PolInSAR coherence. These two decorrelations are directly related to the vertical structure, temporal stability, and physical properties of the forest scattering medium, and the compensated PolInSAR coherence can be further approximated as a coherence function below based on the two-layer dynamic scattering model [15].
where, φ g is the ground phase, γ vt is the coherence component with coupled effects originating from the volumetric scatterer and the wind-derived temporal changes, and µ → ω is the ground-to-volume amplitude ratio varying with polarization → ω. The temporal changes in the ground layer are neglected due to the short temporal baselines of the UAVSAR data. The volumetric-temporal coherence can be further described as follows: where k z is the vertical wavenumber provided by where λ is the radar wavelength, ∆θ is the difference between the repeat-pass incidence angles, B ⊥ is the perpendicular baseline, and R is the slant range. By substituting Equations (1)-(4) into (7), the volumetric-temporal coherences of the physical model with LVM, LVA, QVM, and QVA can be described by Equations (9)- (12).
The volumetric-temporal coherences in Equations (9)- (12) are plotted in the complex plane using the same parameters in Figure 4 (the volume height from 0 to 25 m, center frequency of 1.26 GHz, and the vertical wavenumber of 0.2 for an incidence angle of 40 • ). Each blue segment in Figure 5 corresponds to an RMoG volume with the LVA and LVM of a given height and mean extinction from 0 to 1 dB/m, where the volumes with lower extinctions have smaller radial distances and thus lower coherences. Compared with the blue segments, the introduction of QVA and QVM have little impact on the volumetric-temporal coherence for lower heights and extinctions, as presented in Figure 5b-d. However, larger differences are observed in volumes with higher forest heights and extinctions. Figure 5b illustrates that introducing QVM accelerates the loss of coherence magnitude with increasing height, especially in volumes with higher extinction. As presented in Figure 5c, introducing QVA into the volume leads to a smaller interferometric phase at the same coherence magnitude level as the RMoG volume, especially in the condition of high extinction. Volumetric-temporal coherences derived from the volumes with both QVA and QVM are shown in Figure 5d, where smaller coherence magnitude and interferometric phase are observed at the same time compared with the RMoG volumes.
of 40°). Each blue segment in Figure 5 corresponds to an RMoG volume with the LVA and LVM of a given height and mean extinction from 0 to 1 dB/m, where the volumes with lower extinctions have smaller radial distances and thus lower coherences. Compared with the blue segments, the introduction of QVA and QVM have little impact on the volumetric-temporal coherence for lower heights and extinctions, as presented in Figure 5b-d. However, larger differences are observed in volumes with higher forest heights and extinctions. Figure 5b illustrates that introducing QVM accelerates the loss of coherence magnitude with increasing height, especially in volumes with higher extinction. As presented in Figure 5c, introducing QVA into the volume leads to a smaller interferometric phase at the same coherence magnitude level as the RMoG volume, especially in the condition of high extinction. Volumetric-temporal coherences derived from the volumes with both QVA and QVM are shown in Figure 5d, where smaller coherence magnitude and interferometric phase are observed at the same time compared with the RMoG volumes. Each blue segment in Figure 6 corresponds to an RMoG volume with the LVA and LVM of a given height and the random motion coefficient ranging from 0 to 0.1, where volumes with stronger motions have smaller radial distances. Similar behaviors of the Each blue segment in Figure 6 corresponds to an RMoG volume with the LVA and LVM of a given height and the random motion coefficient ranging from 0 to 0.1, where volumes with stronger motions have smaller radial distances. Similar behaviors of the volumetrictemporal coherences can be observed in Figure 6 as in Figure 5, such as the introduction of QVM has a stronger impact on the coherence magnitude than the interferometric phase. As shown in Figure 6b, coherence magnitudes drop faster with increasing height. Likewise, introducing QVA into the volume leads to a smaller interferometric phase for the same coherence magnitude level. Volumetric-temporal coherences derived from the synergy of QVA and QVM also have smaller coherence magnitudes and interferometric phases at the same time compared with the RMoG volumes, as presented in Figure 6d.
introduction of QVM has a stronger impact on the coherence magnitude than the interferometric phase. As shown in Figure 6b, coherence magnitudes drop faster with increasing height. Likewise, introducing QVA into the volume leads to a smaller interferometric phase for the same coherence magnitude level. Volumetric-temporal coherences derived from the synergy of QVA and QVM also have smaller coherence magnitudes and interferometric phases at the same time compared with the RMoG volumes, as presented in Figure 6d.

Forest Height Inversion
The PolInSAR coherence of a two-layer dynamic scattering model is fully determined by five parameters (ground phase, forest height, random motion factor, groundto-volume amplitude ratio, and the scattering attenuation coefficient). Three of them (scattering attenuation coefficient, random motion factor, and forest height) are accounted for in the polarization-independent volumetric-temporal coherence term equivalent to the PolInSAR coherence in the absence of ground backscatter contributions with the ground phase removed. Equation (6) indicates a linear signature of coherences predicted by the proposed model on the unit circle, which can be represented by a straight line varying with the ground-to-volume ratio in the complex plane. This fitted line intersects the unit circle at two points, and one of them is the correct ground solution. We select the correct one by assuming that the height of the highest observed phase center is less

Forest Height Inversion
The PolInSAR coherence of a two-layer dynamic scattering model is fully determined by five parameters (ground phase, forest height, random motion factor, ground-to-volume amplitude ratio, and the scattering attenuation coefficient). Three of them (scattering attenuation coefficient, random motion factor, and forest height) are accounted for in the polarization-independent volumetric-temporal coherence term equivalent to the PolInSAR coherence in the absence of ground backscatter contributions with the ground phase removed. Equation (6) indicates a linear signature of coherences predicted by the proposed model on the unit circle, which can be represented by a straight line varying with the ground-to-volume ratio in the complex plane. This fitted line intersects the unit circle at two points, and one of them is the correct ground solution. We select the correct one by assuming that the height of the highest observed phase center is less than π/k z m above the ground, as indicated in [12,21]. After that, exp jφ g γ vt can be determined by the furthest coherence point from the estimated ground solution, and the volumetric-temporal coherenceγ vt can be estimated by removing the ground phase [22].
However, a volumetric-temporal coherence estimation can only provide magnitude and phase information with two DOFs and causes underdetermination for the inversion of three unknown parameters. Therefore, an extra pass of the UAVSAR data is acquired to form dual-baseline repeat-pass interferometry. Each interferometry provides a volumetrictemporal coherence estimation through the above-mentioned method. Since repeat-pass UAVSAR data are acquired at slightly different altitudes and incidence angles within hours, we assume a stable attenuation property (a constant attenuation coefficient) in the forest canopy within different passes. In this condition, one additional interferometry can provide volumetric-temporal coherence estimations with two more DOFs while only introducing one additional unknown motion factor. Then, the forest height can be estimated by solving nonlinear equations using iterative optimization to minimize the least square distance between the predicted and observed volumetric-temporal coherences.
where, f i is the volumetric-temporal coherence function. Equations (9)-(12) provide four descriptions of the volumetric-temporal coherence for the forest height inversion. All of them can be employed in (13) as f i and each one of them corresponds to a certain attenuation and random motion profiles in the two-layer dynamic scattering model. In the iterative optimization of each model, a minimized square distance can be obtained regarding the minimum difference between the predicted and observed volumetric-temporal coherences. To leverage various scattering attenuation and random motion profiles, we employed a pixel-wise optimization strategy by choosing forest height solutions derived from the coherence function f i with the smallest minimized square distance in each pixel, which is as follows: This pixel-wise optimization strategy can be interpreted as selecting the best-fit scattering model to depict trees in each pixel, which is much more reasonable than a global scattering model because radar images are mostly with medium resolutions, and thereby, the trees in each pixel may exhibit various behaviors.
These forest height results are quantitatively evaluated by computing the bias and Root Mean Square Error (RMSE) versus the LVIS RH100 trees height, given by:

Scattering Attenuation Fitting
By specifying the permittivities for the forest elements and integrating the amplitudes based on a forward scattering model, attenuation grids can be generated by Pol-SARproSim+ [35]. The effective attenuations derived from the attenuation grid are used in the following inversion analysis. Since these attenuation grids provide estimates of the attenuation at the radar spatial resolution with a one-meter vertical resolution, it maintains the heterogeneity of the tree crown distribution. To assess the applicability of LVA and QVA to a particular tree species, the forest stands and scattering attenuation maps of the 14 deciduous and four coniferous species listed in Table 1 are simulated by tree models embedded in PolSARProSim+. The average attenuation profiles derived from these output attenuation maps are plotted as blue curves in Figures 7 and 8. The red and yellow curves are the best-fit LVA and QVA descriptions compared with the associated blue curve, and the RMSEs of the fitting results are presented in the plots in Figures 7 and 8. LVA and QVA to a particular tree species, the forest stands and scattering attenuation maps of the 14 deciduous and four coniferous species listed in Table 1 are simulated by tree models embedded in PolSARProSim+. The average attenuation profiles derived from these output attenuation maps are plotted as blue curves in Figures 7 and 8. The red and yellow curves are the best-fit LVA and QVA descriptions compared with the associated blue curve, and the RMSEs of the fitting results are presented in the plots in Figures 7 and 8.   Figure 7. Vertical scattering attenuation fitting results of the deciduous forest stands (see Table 1 for an explanation of the acronyms).  Table 1 for an explanation of the acronyms). Remote Sens. 2022, 14, x FOR PEER REVIEW 13 of 20 Figure 8. Vertical scattering attenuation fitting results of the coniferous forest stands (see Table 1 for an explanation of the acronyms).
The attenuation description with a smaller RMSE is chosen as the best-fit attenuation model for the forest stands of a certain species, as listed in Table 2, indicating that LVA is more suitable to describe coniferous species while a mixture of LVA and QVA is needed for different deciduous species.

Boreal Ecosystem Research and Monitoring Sites
Dual-baseline PolInSAR data are selected from multi-baseline UAVSAR collections over BERMS. The vertical baselines between the primary and two auxiliary images are 40 m, and the temporal baselines are 25 min, with an average absolute vertical wavenumber of 0.09 and a / of 35 m. Since the maximum forest height in BERMS is about 33 m, the selection of longer vertical baselines will result in a smaller ambiguity height than the canopy height and hence an invalid assumption during the ground solution selection. The forest height results derived from the UAVSAR data using the reference model with various attenuation and motion properties are presented in Figure 9, where areas with low backscatter (less than −15 dB for an incidence angle from 0-35°, less than −21 dB for an incidence angle from 35-45°, less than −24 dB for an incidence angle from 45-55°, and less than −28 dB for an incidence angle larger than 55°) are masked out [38] Table 1 for an explanation of the acronyms).
The attenuation description with a smaller RMSE is chosen as the best-fit attenuation model for the forest stands of a certain species, as listed in Table 2, indicating that LVA is more suitable to describe coniferous species while a mixture of LVA and QVA is needed for different deciduous species. Dual-baseline PolInSAR data are selected from multi-baseline UAVSAR collections over BERMS. The vertical baselines between the primary and two auxiliary images are 40 m, and the temporal baselines are 25 min, with an average absolute vertical wavenumber of 0.09 and a π/k z of 35 m. Since the maximum forest height in BERMS is about 33 m, the selection of longer vertical baselines will result in a smaller ambiguity height than the canopy height and hence an invalid assumption during the ground solution selection. The forest height results derived from the UAVSAR data using the reference model with various attenuation and motion properties are presented in Figure 9, where areas with low backscatter (less than −15 dB for an incidence angle from 0-35 • , less than −21 dB for an incidence angle from 35-45 • , less than −24 dB for an incidence angle from 45-55 • , and less than −28 dB for an incidence angle larger than 55 • ) are masked out [38] and a 3 × 3 moving average is performed to smooth the forest height maps. The related PolInSARderived heights are plotted versus the LVIS RH100 height for a quantitative comparison, as shown in Figure 10, and the associated RMSE, BIAS, and R 2 metrics are given in Table 3. As indicated in Section 2.2, the spatial resolutions of the retrieved forest height map are 4.8 m in azimuth and 3.2 m in the slant range, but the RH100 metrics derived from the LVIS full-waveform LiDAR measurements are much sparser. Thereby, the down-sampling process is conducted over the retrieved forest height for a quantitative comparison, where the UAVSAR-derived height samples are uniformly selected from the whole forest height map with a separation of 10 pixels in both the slant range and azimuth directions, and the pixels corresponding to the null height values in the RH100 map are not included in the calculation of any accuracy metrics.
Remote Sens. 2022, 14, x FOR PEER REVIEW 14 of 20 and a 3 × 3 moving average is performed to smooth the forest height maps. The related PolInSAR-derived heights are plotted versus the LVIS RH100 height for a quantitative comparison, as shown in Figure 10, and the associated RMSE, BIAS, and R 2 metrics are given in Table 3. As indicated in Section 2.2, the spatial resolutions of the retrieved forest height map are 4.8 m in azimuth and 3.2 m in the slant range, but the RH100 metrics derived from the LVIS full-waveform LiDAR measurements are much sparser. Thereby, the down-sampling process is conducted over the retrieved forest height for a quantitative comparison, where the UAVSAR-derived height samples are uniformly selected from the whole forest height map with a separation of 10 pixels in both the slant range and azimuth directions, and the pixels corresponding to the null height values in the RH100 map are not included in the calculation of any accuracy metrics.   Dual-baseline PolInSAR data are also selected from the multi-baseline collection over the PNP site, as described previously in Section 2.2. The maximum forest height in the PNP site is much higher than the BERMS (~65 m) and thus requests for a smaller vertical baseline with a larger ambiguity height. The vertical baselines between the primary and two auxiliary images are both 20 m, with temporal baselines of 25 min and 50 min, respectively. The related average absolute vertical wavenumber is about 0.05 with / of ~63 m. Forest height results over the PNP site based on the reference model in Section 3.1 are shown in Figure 11. Similar masking and smoothing operations are conducted over the PNP site, as was conducted for BERMS. The density plots between the estimated forest heights and the LVIS RH100 height are presented in Figure 12, along with the associated RMSE, BIAS, and R 2 metrics given in Table 4.

Pongara National Park
Dual-baseline PolInSAR data are also selected from the multi-baseline collection over the PNP site, as described previously in Section 2.2. The maximum forest height in the PNP site is much higher than the BERMS (~65 m) and thus requests for a smaller vertical baseline with a larger ambiguity height. The vertical baselines between the primary and two auxiliary images are both 20 m, with temporal baselines of 25 min and 50 min, respectively. The related average absolute vertical wavenumber is about 0.05 with π/k z of~63 m. Forest height results over the PNP site based on the reference model in Section 3.1 are shown in Figure 11. Similar masking and smoothing operations are conducted over the PNP site, as was conducted for BERMS. The density plots between the estimated forest heights and the LVIS RH100 height are presented in Figure 12, along with the associated RMSE, BIAS, and R 2 metrics given in Table 4.  Dual-baseline PolInSAR data are also selected from the multi-baseline collection over the PNP site, as described previously in Section 2.2. The maximum forest height in the PNP site is much higher than the BERMS (~65 m) and thus requests for a smaller vertical baseline with a larger ambiguity height. The vertical baselines between the primary and two auxiliary images are both 20 m, with temporal baselines of 25 min and 50 min, respectively. The related average absolute vertical wavenumber is about 0.05 with / of ~63 m. Forest height results over the PNP site based on the reference model in Section 3.1 are shown in Figure 11. Similar masking and smoothing operations are conducted over the PNP site, as was conducted for BERMS. The density plots between the estimated forest heights and the LVIS RH100 height are presented in Figure 12, along with the associated RMSE, BIAS, and R 2 metrics given in Table 4.

Boreal Ecosystem Research and Monitoring Sites
Compared with the RH100 metrics, the forest heights derived from the reference model over the BERMS site are generally overestimated, as shown in Figures 9 and 10. A part of this error is related to the mismatch of acquisition times between the UAVSAR and LVIS data. However, it can be concluded from the difference between the RH100

Boreal Ecosystem Research and Monitoring Sites
Compared with the RH100 metrics, the forest heights derived from the reference model over the BERMS site are generally overestimated, as shown in Figures 9 and 10. A part of this error is related to the mismatch of acquisition times between the UAVSAR and LVIS data. However, it can be concluded from the difference between the RH100

Boreal Ecosystem Research and Monitoring Sites
Compared with the RH100 metrics, the forest heights derived from the reference model over the BERMS site are generally overestimated, as shown in Figures 9 and 10. A part of this error is related to the mismatch of acquisition times between the UAVSAR and LVIS data. However, it can be concluded from the difference between the RH100 metrics of 2017 and 2019 (RMSE: 0.58 m) that this temporal error only contributes a small part to the total. The other part of the error is associated with the mismatch between the forest scenario and the reference model used, where the ground scattering contribution is neglected during the height inversion. This neglection of ground scattering contributes more to the total error compared with the temporal mismatch because the overestimation is much more prominent in smaller trees than the taller ones. QVM in the reference model significantly reduces the overestimation and exhibits a better performance in either the LVA or QVA conditions, as illustrated in Figure 10b,d. A 1.75 m reduction in the RMSE and a 2.38 m reduction in the bias are observed in LVA + QVM results, and a 2.13 m reduction in the RMSE and a 2.61 m reduction in bias are achieved by the QVA + QVM result. This suggests that a Gaussian distribution with zero mean and vertically quadratic variance better fits the dynamic motion properties of forests in BERMS at the UAVSAR data acquisition time. In addition to that, the reference model with LVA surpasses QVA for forest height inversion in either LVM or QVM conditions because the heterogeneous attenuation aggravates the overestimation of smaller trees, as indicated in Figure 10c,d. A 1.25 m reduction in RMSE and a 1.31 m reduction in bias are achieved by the LVA + LVM result, and a 0.87 m reduction in RMSE and a 1.08 m reduction in bias are observed in the LVA + QVM result. The advantages of selecting a proper scattering attenuation and random motion structure for the forest height inversion can be observed when using the pixel-wise optimization strategy, which leverages LVA, LVM, QVA, and QVM on a pixel basis and, therefore, achieves better performance than others that do not vary spatially. As presented in Figures 9f and 10e, a 0.35 m reduction in the RMSE and a 0.60 m reduction in bias are observed compared with the best-performing single model (LVA + QVM).

Pongara National Park
Compared with the RH100 metrics, the forest heights derived from the reference model over the PNP site are generally underestimated, as presented in Figures 11 and 12, particularly in areas with taller trees. One possible explanation for this underestimation is the overcompensation of temporal decorrelation because the underestimation is more prominent in results derived from the reference model with QVM. Contrary to boreal forests, tropical rainforests experience very light winds and, therefore, fewer temporal changes occur in the ground layer. The LVM description in the reference model significantly reduces the underestimation and exhibits a better forest height estimation performance in either LVA or QVA conditions, as presented in Figure 12a,c. A 3.12 m reduction in RMSE and a 4.55 m increase in bias are achieved by the LVA + LVM result, whereas a 1.26 m reduction in the RMSE and a 4.68 m increase in bias are observed in the QVA + LVM result, indicating that a Gaussian distribution with zero mean and vertically linear variance better fits the dynamic motion properties of tropical forests in PNP. Moreover, the reference model with QVA performs better on the forest height inversion than the LVA ones in either LVM or QVM conditions by reducing the underestimation of taller trees, as shown in Figure 12c,d. A 0.88 m reduction in RMSE and a 3.17 m increase in the bias are obtained by the QVA + LVM result, and a 2.74 m reduction in the RMSE and a 3.04 m increase in the bias are achieved by the QVA + QVM result. Figures 11f and 12e also give evidence to support the superiority of the pixel-wise optimization, where a 0.35 m reduction in the RMSE and a 0.02 m reduction in the bias are achieved compared with the best-performing single model (QVA + LVM).

Model Analysis
It can be concluded from above that the model with a homogeneous attenuation profile (LVA) and a heterogeneous motion profile (QVM) has a better forest height inversion performance in boreal forest such as BERMS. Meanwhile, the model with a heterogenous attenuation profile (QVA) and a homogeneous motion profile (LVM) achieve better results in tropical forests such as PNP. These conclusions are directly related to wind conditions and tree species in both forest sites.
On the one hand, according to the historical weather and climate data achieved by the government of Canada [39], the wind speed at the data acquisition time is around 12 km/h over the BERMS site, making leaves move and small branches sway. The wind is much lighter over PNP, which is less than 6 km/h with reference to the global historical weather and climate data [40]. Compared with PNP, the stronger wind in BERMS makes different parts of trees (leaves, stem, and small branches) participate in the dynamic process and thereby related temporal decorrelations are more sensitive to the intrinsic forest structure, and a heterogenous description of the random motion achieved a better performance in the forest height inversion.
The attenuation property, on the other hand, is more related to tree species in the forest sites. This has been indicated in the scattering attenuation fitting results derived from the PolSARproSim+ simulated data (see Section 4.1), where the model with LVA is more suitable to describe the coniferous species, and the model with LVA or QVA is more suitable for the deciduous species. Since the BERMS forest mainly consists of coniferous species such as jack pine, black spruce, and aspen, the reference model with LVA has a better performance in forest height inversion in this area. Similarly, the PNP forest is mostly covered by deciduous red mangroves (RHMA), and thus, the reference model with QVA achieves better results in this tropical area.

Conclusions
In this paper, we demonstrated the potential of dual-baseline repeat-pass airborne PolInSAR for estimating forest height at the landscape scale. A physical model with LVA, QVA, LVM, and QVM profiles in the volume layer is proposed based on the RMoG model to investigate the impacts of homogeneous or heterogeneous attenuation and random motion properties on the performance of forest height inversion. The experiments were carried out over a boreal and tropical forest site to show the efficacy of different volumetric profiles over various forest types. The related forest height results were compared with LVIS RH100 references for the quantitative evaluation. Volumes with LVA and QVM achieved the best performance (RMSE of 3.56 m and bias of 2.05 m) over the boreal forest dominated by coniferous tree species, while volumes with QVA and LVM exhibited the best performance (RMSE of 6.83 m and bias of 0.43 m) in the tropical forests dominated by deciduous red mangroves. This is consistent with the scattering attenuation fitting results derived from the PolSARproSim+ simulated data. The forest height results generated by the pixel-wise optimization strategy surpass the best-performing single models in both the boreal (RMSE of 3.21 m and bias of 1.45 m) and tropical (RMSE of 6.48 m, bias of 0.41 m) forests, which indicates its superiority in leveraging the advantages of different volumetric profiles.