A Forest Height Joint Inversion Method Using Multibaseline PolInSAR Data

Estimating vegetation height from polarimetric interferometric synthetic aperture radar (PolInSAR) data using the random volume over ground (RVoG) model has long been used. Most of these methods propose models and apply them to real airborne data to demonstrate their potential. The single-baseline PolInSAR forest height estimation based on the RVoG model lacks sufficient observation information. For this reason, multibaseline data are introduced to address this. This letter fits the relationship of model parameters in multibaseline observation scenarios and focuses the forest height inversion on the calculation of pure volume decorrelation. Subsequently, a multibaseline forest height joint inversion method based on the least-squares principle is adopted. Finally, we use airborne PolInSAR data from the Lope and Mondah sites collected by uninhabited aerial vehicle synthetic aperture radar (UAVSAR) and F airborne synthetic aperture radar (F-SAR) systems during AfriSAR 2016 to verify the proposed method. The experimental results show that the accuracy of the proposed method (Lope: root mean square error (RMSE) = 5.8 m, Mondah: RMSE = 5.12 m) is 38.1% and 34.53% higher than the coherence separation product (Lope: RMSE = 9.37 m, Mondah: RMSE = 7.82 m).


I. INTRODUCTION
F OREST height determination is of great importance for the management of forest resources and the investigation of the global carbon cycle [1], [3], [5]. Optical remote sensing, light detection and ranging (LiDAR), and synthetic aperture radar (SAR) have been used for measuring forest height. Polarimetric interferometric synthetic aperture radar (PolInSAR), due to its spatial continuity and large coverage, is regarded as a potential tool to extract global forest height [1], [2], [5], [6].
To extract forest height from PolInSAR data, it is necessary to find a scattering model that can describe the scattering process of SAR signals penetrating the forest layer. To achieve this goal, the random volume over ground (RVoG) model that combines PolInSAR observations with forest biophysical parameters has been proposed [1], [2], [3]. Based  RVoG model, Papathanassiou and Cloude [3] and Cloude and Papathanassiou [4] employed a nonlinear iteration and three-stage algorithm to retrieve forest height from singlebaseline PolInSAR data. Yet, in the absence of sufficient observations, both these algorithms need to assume that there is a polarization channel whose ground-to-volume ratio (GVR) is less than −10 dB [4]. This restriction cannot be met in many cases. Although a modified three-stage method, the fixed extinction method [5], has been used to solve this issue, it is difficult to obtain a reliable extinction coefficient in practice.
In addition, spatially inconsistent forest height distributions usually require PolInSAR data with different interferometric baselines to meet height sensitivity requirements. The above two problems motivate us to use multibaseline data to retrieve forest height. Multibaseline algorithms can be broadly classified into three categories. The first category considers the effect of vertical wavenumber on the forest height estimation to select an interferometric baseline [6], [7], [8], [17], [19]. However, in fact, only single-baseline PolInSAR data is used for every resolution cell and null GVR assumption is still needed. To solve this problem, based on the dual-baseline inversion algorithm [9], Liao et al. [13] proposed a multibaseline PolInSAR inversion algorithm. However, this method cannot consider the influence of the vertical wavenumber on forest height estimation. The third category is the comprehensive utilization of multibaseline PolInSAR data, such as multibaseline nonlinear iteration, which can avoid the null GVR assumption and consider the impact of vertical wavenumber on the forest height estimation [14], [21]. However, there are still two remaining problems that should be considered: 1) the formula of the RVoG model is complex, reliable initial values for all the unknown parameters are needed to avoid local-optimal problem, which is not easy in practice, especially for forest height and extinction [21] and 2) the existing multibaseline algorithm focuses on scaling the effect of vertical wavenumber on the forest height inversion. However, for PolInSAR observations combing with the RVoG model, the main factor that determines the forest height accuracy is the errors in PolInSAR observations [20], [22].
To consider the above problems, a new forest height estimation method based on multibaseline PolInSAR data is proposed. To reduce the complexity of the RVoG inversion, our method can be divided into two steps. First, pure volume coherence estimation: when using the RVoG model, the multibaseline PolInSAR data are combined by considering the linear relationship among phases of volume decorrelations corresponding to different baselines, and the Cramer-Rao lower bound is used to describe the errors in PolInSAR observations. The volume decorrelation is then estimated by 1558-0571 © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
using the least-squares-based method. Second, forest height estimation: The forest height and extinction are retrieved from the estimated pure volume decorrelation by a 2-D lookup table [4], [13]. This letter is structured as follows. Section II introduces the linear constraint between the volume decorrelation of the multibaseline RVoG model and the forest height inversion method. Section III mainly derives the forest height results and discussion. Finally, the conclusions are drawn in Section IV.

A. Relationship Among Different Volume Decorrelations
The RVoG model has been widely used as a coherent scattering model for forest parameter retrieval. It correlates PolInSAR observations with forest parameters as follows [1], [3], [4]: where w is the polarization vector, ϕ 0 is the ground phase, and μ(w) represents the GVR, which changes with the polarization channel. γ v is the volume-only decorrelation caused by the vegetation layer, which can be written as [1], [3] γ where h v is forest height and σ represents the wave extinction coefficient of the vegetation layer. θ is the incident angle. k z is the vertical wavenumber of the interferometer. For multibaseline interferograms, the constraint relationship among volume decorrelations should be considered. According to formula (2), there is an obvious linear relationship in the change in the pure volume decorrelation phase α with k z [16]. The pure volume decorrelation phase relations for different interferometric baselines are fit as The superscript i = 1, 2, . . . , n represents the baseline number. Formula (3) is simplified, and only the two most basic variables α 1 and e are retained to represent other phase values α i It should be noted that this fitting is approximate. To determine whether this relationship can be used in the inversion of forest height, we use formula (4) to fit phase α i fit and then compare it with the phase α i true of γ v in formula (2). The fitting deviation can be expressed as The results show that when the extinction is 0.1 dB/m, the fitting line of the volume decorrelation phases has the maximum deviation, and the maximum error is 0.71 m when the forest height is 30 m. For the other extinction in formula (2), the fitting error is less than 0.3 m. This amount of deviation is acceptable in forest height inversion.

B. Calculation of Volume Decorrelations
For the RVoG model in (1), the complex interferometric coherence can be expressed as γ (ω) = R(ω) + i I (ω), and the ground component can be written as e iϕ 0 = cos ϕ 0 + i sin ϕ 0 . The volume component γ v can be written as A(cos α+i sin α), where A and α are the magnitude |γ v | and phase angle(γ v ) of the volume decorrelation, respectively. Then, the model in formula (1) can be rewritten as where superscript n represents the baseline number. For the same target observed at multibaseline, μ(w) with the same polarization channel from different baselines can be considered to be the same [18]. In contrast, the volume decorrelation coefficient γ v and the terrain phase ϕ 0 do not change with the polarization but vary with different baseline observations. In the multibaseline calculation, we add formulae (4)-(6) and then solve the multibaseline RVoG model using the principle of the least squares [14], [15]. Simultaneously, a Taylor series expansion is performed on the real and imaginary parts, and only the zero-order and first-order terms are retained For formula (7), this can be summarized in the form of the general formula in the multibaseline inversion v 2nm×1 = B 2nm×(2n+m+2) ·x (2n+m+2)×1 + l 2nm×1 (8) where m and n represent the number of polarization channels and baselines, respectively. For the observation value v, each polarization complex coherence coefficient corresponds to two observation values. Therefore, the number of observed values is 2 nm. The parameterx to be determined includes the ground phase ϕ 0 , the amplitude A, the GVR μ(w), and the α 1 and e, for a total of 2n + m + 2 parameters.
To keep the observation data as independent as possible, this letter uses the maximum phase separation optimization algorithm to obtain the two complex coherence coefficients with the largest phase difference [11], [12]. In this way, the number of polarization channels is m = 2. The baseline data with a large difference from the model can be eliminated (|l| > ε). The value of ε is set to 1. According to formula (8) and the least-squares principle, v T Pv = min, the value of the parameterX to be determined is [14], [15] x = B T PB −1 BPl;X = X 0 +x.
For the initial value of the parameter X 0 in the solution process, the ground phase approximation is estimated by the fitting straight line of the polarization complex coherence coefficient on the complex plane circle [2], [3], [4]. Then, we assume that the GVR of the polarization channel PDHigh μ h = 0 as an approximation. The volume decorrelation and μ l are further calculated according to formula (1). These approximations are then substituted into the multibaseline RVoG model [formula (7)] to obtain the coefficient matrices B and l. The adjustment valuex and the true value (unbiased value)X can then be calculated according to formula (9). P represents the weight function that adjusts the contribution of different interference baselines, which is described in detail in Section II-C.

C. Weight Function
Different interferograms may produce different levels of forest height deviation [3]. To reasonably solve the least-square function, accurate prior weight is essential. In this letter, the standard deviation of the coherence amplitude (Cramer-Rao bound) is used to determine the weight according to the absolute value of the observed complex coherence coefficient [15]. The Cramer-Rao lower bound on the variance of coherence magnitude can be approximately related to the coherence of different baselines and polarizations as [11] where L is the number of independent samples used to estimate the coherence. We assume that observations between different polarizations and different baselines are independent when using multibaseline PolInSAR data. The weights of the real part R(ω) and the imaginary part I (ω) of the same baseline are the same and independent. The variance of the polarization with the smallest error is selected as the unit weight variance. The higher the interferometric coherence is, the smaller the variance and the larger the weight P(w i ). The corresponding expression is

D. Forest Height Estimation
Through this adjustment, the solution of the parameters is closer to the true value than the approximate value given by the three-stage algorithm. Due to the constraints of the weight function, the improvement of the parameters in the quality of the data with smaller errors is greater. According to the principle of least squares, a matrix Qxx = (B T P B) −1 can be used to judge the improvement. The smaller the Q value is, the smaller the difference between the solution and the true value [15]. We select the volume coherence γ i v corresponding to the smallest difference value to calculate the forest height. Putting this into practice, we find that the solution of the volume coherence magnitude A is the most robust, which is due to its small range of values. The baseline number i corresponds to the smallest Q A value can then be selected from among all baselines Next, according to formula (2), the inversion problem of forest height can be simplified to a 2-D nonlinear problem. The solution can be obtained using a 2-D lookup table by selecting the forest height and extinction corresponding to the   smallest difference in the volume coherence coefficient [4]. The process of the multibaseline joint calculating forest height based on the least-squares principle is shown in Fig. 1.

III. EXPERIMENTAL RESULTS AND DISCUSSION
In order to test the proposed method, two test sites, Lope and Mondah, were selected. Two test sites are covered with typical African tropical forest, with heights ranging from 5 to 55 m. Airborne PolInSAR data were collected by the AfriSAR 2016 campaign. The baseline parameters of different interferometric pairs are listed in Table I. The Land, Vegetation, and Ice Sensor (LVIS) lidar forest height product (RH100) was obtained by NASA during this campaign. The resolution of LVIS forest height product is 25 × 25 m. For the accuracy validation, we resampled (linear interpolation) it to the same resolution as the multilooked SAR image (see Fig. 2).
The L-band PolInSAR data over the Lope test site was acquired by the NASA uninhabited aerial vehicle synthetic aperture radar (UAVSAR) system. The original polarimetric SAR image has a resolution of 0.6 m in azimuth and 1.67 m in the slant range. We performed multilooking using a rectangular window with a size of 16 pixels in azimuth and 4 pixels in slant range, producing multilooked images with a pixel spacing of 9.6 m in azimuth and 6.66 m in slant range. The main characteristics of multilooking are resample and mean window filter. Six interferometric pairs were selected for testing the proposed method. For the coherent separation algorithm in Lope, when the extinction coefficient is 0.1 dB/m, the overall inversion accuracy is optimal, and the center of the scatter plot passes through the line y = x. The coherence separation product with a fixed extinction coefficient σ = 0.1 dB/m was also used to estimate the forest height [6], [8], so that we can make a comparison with the proposed method. The final results are shown in Fig. 3. It is clear that, with respect to the LVIS RH100 forest height product, the proposed multibaseline joint solution method gives a more accurate forest height estimation result than the coherent separation product in detail. This confirms that it is necessary to consider different GVR and extinction coefficients for each observation cell.
For the more intuitive accuracy assessment, we validate the inversion results with LVIS lidar data, as shown in Fig. 4(a) and (b). A height less than 3 m was considered grassland and masked. For the PolInSAR and LVIS forest height image, we divide it into blocks of 40 × 40 windows. The average forest height in 40 × 40 window (260 × 260 m) was used for point-to-point statistics, and a total of 10 734 forest statistical blocks were counted. For the forest height obtained by the coherence separation product, its root mean square error (RMSE) was 9.37 m, and the correlation coefficient (R 2 ) was 0.74. The RMSE of the proposed method is 5.8 m with an R 2 of 0.89. This can be noted that PolInSAR forest heights are overestimated in low vegetation and underestimated in high vegetation compared to LVIS forest height. This can be explained: For the L-and P-band PolInSAR data, it is difficult to select a polarization channel without recording the ground scattering contribution, because of its significant penetration depth. In order to solve this problem, in the optimal coherence separation method, the extinction of the RVoG model is set to a known constant which is identical for the whole image. This results in the line y = x passing through the center of the scatter plot when the best accuracy is achieved. In the joint solution method, multibaseline data is introduced to solve the GVR and extinction. More interferometric baselines, however, mean that some are not suitable for estimating different forest heights. In low forests, the introduction of longer baseline data leads to an overestimation; in high forests, some shorter baselines lead to an underestimation of it [16], [20]. Since the error of the observed data is reduced by the least-squares method, the magnitude of this forest height estimation bias is much smaller than that of the coherent separation product. This confirms that in the absence of observation errors, all interferometric baselines can correspond to an accurate forest height, and the influence of vertical wavenumber on inversion performance can be ignored.
The P-band PolInSAR data over Mondah was collected by the German Aerospace Center (DLR)'s F airborne SAR (F-SAR) system. The resolution of the PolSAR image is 0.81 m in azimuth and 1.2 m in the slant range. Multilooking processing was performed for every 4 pixels in azimuth and every 2 pixels in the slant range, producing multilooked images with a pixel spacing of 3.24 m in azimuth and 2.4 m in the slant range. For P-band data in Mondah, when the extinction coefficient is 0.4 dB/m, it has the best overall inversion accuracy. The coherence separation product with a fixed extinction coefficient σ = 0.4 dB/m was used for the comparison of the methods. The final PolInSAR forest height maps are shown in Fig. 5. We used the average forest height in a window of size 40 × 40 (96 × 96 m) pixels for the point-to-point statistics, and a total of 6659 forest statistical blocks were used for the statistics. The statistical results are shown in Fig. 4(c) and (d), where the RMSE of the coherence separation product is 7.82 m, and the R 2 is 0.76. The RMSE of the proposed method is 5.12 m with R 2 of 0.91.
The baseline distribution used in the experimental verification is relatively uniform because the baseline distribution obtained by the AfriSAR campaign is uniform, which is mainly limited by the form of data. However, the method proposed in this article has nothing to do with the uniform distribution of the baseline. In the Mondah test site, we selected a set of incompletely uniform baseline distributions for reverification. Through the experimental comparisons, it can be seen that the proposed method has stable improvement in both experimental areas. Especially in the low vegetation area, it is more obvious. This is because the multibaseline data makes the solution of GVR more accurate and eliminates the influence of ground scattering. Second, compared with the coherence separation product, the correlation coefficient is greatly improved, which proves that the least squares can eliminate the random error when it is applied to the calculation of the pure volume decorrelation. Third, the proposed joint solution method, under the control of the weight function, gives different weights to different vertical wavenumbers. This can select the observation data with small error, which can avoid the influence of abnormal observation data on forest height estimation. The method of multibaseline data joint inversion improved the overall deviation in forest height product and reduced the requirements of the RVoG model for forest height inversion conditions.

IV. CONCLUSION
In this letter, we have established a constraint relationship among different volume decorrelation coefficients and proposed a robust joint solution method for the multibaseline RVoG model to estimate forest height. This method divides the multibaseline inversion problem into two steps: pure volume decorrelation coefficients solution and forest height estimation. Using this method, the accuracy of the solutions is improved for different baseline data. To demonstrate the quality of the multibaseline joint adjustment method, AfriSAR 2016 airborne L-band and P-band multibaseline data are used for experimental verification. The results show that the proposed method has better inversion performance than the coherence separation product. In future research, more accurate model constraints and temporal decorrelation will be considered to improve the applicability of the inversion method.