Forest Height and Underlying Topography Inversion Using Polarimetric SAR Tomography Based on SKP Decomposition and Maximum Likelihood Estimation

: The key point of forest height and underlying topography inversion using synthetic aperture radar tomography (TomoSAR) depends on the accurate positioning of the phase centers of different scattering mechanisms. The traditional nonparametric spectrum analysis methods (such as beamforming and Capon) have limited vertical resolution and cannot accurately distinguish closely spaced scatterers. In addition, it is very difﬁcult to accurately estimate the ground or canopy heights with single polarimetric SAR images because there is no guarantee that the vertical proﬁle will generate two clear and separate peaks for all resolution cells. A polarimetric TomoSAR method based on SKP (sum of Kronecker products) decomposition and iterative maximum likelihood estimation is proposed in this paper. On the one hand, the iterative maximum likelihood TomoSAR method has a higher vertical resolution than that of the traditional methods. On the other hand, the separation of the canopy scattering mechanism and the ground scattering mechanism is conducive to the positioning of the phase centers. This method was applied to the inversion of forest height and underlying topography in a tropical forest over the TropiSAR2009 test site in Paracou, French Guiana with six passes of polarimetric SAR images. The inversion accuracy of underlying topography of the proposed method was up to 1.489 m and the inversion accuracy of forest height was up to 1.765 m. Compared with the traditional polarimetric beamforming and polarimetric capon methods, the proposed method greatly improved the inversion accuracy of forest height and underlying topography.


Introduction
Forests are important parts of the global ecosystem and play a vital role in the global carbon and oxygen cycle. Long-wavelength synthetic aperture radar has strong penetrating power and can be used to analyze the vertical structure of the forests. Polarization interference synthetic aperture radar (PolInSAR) technology is one of the most important tools for forest vertical structure inversion [1][2][3]. Related studies have proposed a three-stage algorithm [4] and nonlinear iterative algorithm [5] based on the random volume over ground (RVoG) model [4]. Further, they have been successfully applied to the inversion of forest height. However, the performances of PolInSAR methods rely on reasonable model assumptions for the forest, and the calculation of related parameters is complicated.
Synthetic aperture radar tomography (TomoSAR) technology obtains high vertical resolution by synthesizing the aperture in the elevation direction so that it can identify and distinguish scatterers of different heights within the same resolution cell [6]. It has been widely used to reconstruct forest three-dimensional structures and estimate above-ground biomass (AGB) [6][7][8][9][10][11][12][13].
Spectral analysis methods are commonly used in forest structure SAR tomography, including parametric spectral analysis methods (MUSIC) [14] and nonparametric spectral analysis methods (beamforming, Capon) [15]. The main advantages of these methods are their simple algorithm and low computational burden. However, the spectrum analysis methods mentioned above have their own limitations. The MUSIC method needs to know the number of scatterers in a resolution cell, and thus it is not suitable for the estimation of coherent scatterers in forest. The beamforming method has a low vertical resolution and is susceptible to sidelobes. The Capon method has a higher vertical resolution but its radiation resolution is lower [16]. In addition, those spectral analysis methods are seriously affected by the number of SAR images. Compressive sensing (CS) [17][18][19][20][21] can break through the limitation of the number of images and realize high-resolution tomography. However, this method is complicated to operate and has a heavy computational burden [16]. In addition, the performance of the CS method is susceptible to user parameters and thus becomes unstable [12]. Some methods based on statistical regularization and maximum likelihood estimation were proposed for enhanced TomoSAR imaging [22,23]. However, their potential of forest height and underlying topography inversion has not yet been verified, especially in dense tropical forests.
Polarimetric synthetic aperture radar (PolSAR) is sensitive to the shape, direction, and dielectric properties of the scatterers, and it has been widely used in SAR tomography of forests [24][25][26][27]. In [24,25], the authors proposed classic polarization spectrum estimation methods and a high-resolution polarization SAR tomography method, which have been successfully used in under-foliage object imaging. For tomography problems with a small number of measurements, some polarimetric estimators were proposed under the frame of CS and sparsity-based reconstruction [26,27]. These works mainly focused on the reconstruction of the forest profile and the analysis of the scattering mechanism.
The important premise of forest height and underlying topography inversion using SAR tomography is the distinction and accurate positioning of the phase centers of the canopy scattering and ground scattering. However, due to the complexity of the forest structure, the denseness of forest, and the low resolution of the TomoSAR algorithms, there is no guarantee that the vertical profile generates two clear and separate peaks for all resolution cells. Traditional methods have limited resolution; the phase center is not focused enough. Thus, it is difficult to locate the phase centers of different scattering mechanisms. These unfavorable factors can cause errors in the extraction of forest vertical structure parameters.
Some existing studies have shown that obtaining the covariance matrices of different scattering mechanisms is more conducive to the location of the corresponding scattering centers [28][29][30][31]. In order to obtain a more reliable estimation of forest height and underlying topography, a method using polarimetric SAR tomography based on SKP (sum of Kronecker products) decomposition and iterative maximum likelihood estimation (SKP-IMLE) is proposed in this paper. Through SKP decomposition, the forest scattering signal is divided into canopy scattering contribution and ground scattering contribution. Then, the IMLE method is used to estimate the tomograms of different scattering mechanisms. Finally, the corresponding phase center heights are extracted from the tomograms, and the accurate forest height and underlying topography are estimated.
The rest of this paper is organized as follows. Section 2 is devoted to explaining the 3D focus model of TomoSAR and the proposed SKP-IMLE Pol-TomoSAR method. In Section 3, there is an overview of the study area and data set. In Section 4, the effectiveness of the proposed method was evaluated from the real P-band full polarimetric focused SAR images. At the same time, the inversion results of different polarimetric TomoSAR methods were compared. In Section 5, the performance of the proposed method was discussed and analyzed. Finally, conclusions were drawn in Section 6.

TomoSAR Imaging Model
Assuming that the SAR sensor makes N repeated observations on multiple approximately parallel orbits and some pre-processing (such as coregistration, deramping, and phase error calibration [32]) are done, the nth focused complex value for a specific rangeazimuth pixel (r,a) can be established as follows: g(r, x, b n ) = ∆s γ(r, x, s) e jk z (n)s ds (1) where s is the elevation value of the scatterer, and γ(r, x, s) represents the reflectivity function in the elevation direction. k z (n) is the vertical wave number: where b(n) is the vertical baseline, λ is the wavelength and θ is the incident angle. For a multibaseline TomoSAR data set, after D sampling in the elevation direction and considering the noise, (1) can be expressed as the following matrix form: where A = e jk z s d , s d is the height of the scatterer at the dth sampling point, and e is a N-dimensional noise vector.

IMLE TomoSAR Method
Consider that the noise in (3) is the Gaussian white noise with a same variance. The covariance matrix of multibaseline SAR data can be expressed as where (·) H is the conjugate transpose operation. K = diag(k D d=1 ) is the backscatter power matrix, and δ 2 is the noise power and diag(·) converts a vector to a diagonal matrix.
The multibaseline observation data g can be regarded as a complex Gaussian distribution, and its probability density function is defined as Now the maximum likelihood estimation is performed according to the minimum risk Bayesian strategy [23,33]. Ignoring the constant terms, the SAR tomography maximum likelihood estimation problem is equivalent to minimizing the following objective function: Thus, the maximum likelihood estimation of the backscatter power vector k can be regarded as the following minimization problem: According to the solution strategy proposed in [23,33,34], the minimization problem in (7) can be solved when the first derivative of the objective function with respect to k is equal to zero, which leads to where M = diag(k)A H R −1 and (·) diag extracts the diagonal elements of a matrix.
Then, the solution of the IMLE TomoSAR estimator can be organized as follows: whereR is the sample covariance matrix of the multibaseline SAR data. The way to calculate the covariance matrix with the samples in L multilooks is through the following equation: The IMLE method optimizes k through an iterative processing. The initial backscatter power can be estimated from the beamforming method [15]. In addition, the noise variance will be treated as a (diagonal-loading) degree of freedom.

SKP Decomposition
SKP decomposition is an algebraic synthesis method, which can decompose the forest backscattering signal into two parts: the ground scattering contribution (whose phase centers are located on the ground) and the canopy scattering contribution (whose phase centers are located on the canopy) [29,30]. Pol-TomoSAR technology based on SKP decomposition is conducive to the inversion of forest height and underlying topography.
For multibaseline and multipolarization (MBMP) SAR data, the covariance matrix can be expressed as follows: where W denotes the covariance matrix of MPMB data and y = (g T HH g T HV g T VV ) T is the MPMB observation data. E(·) is the operation for mathematical expectation and ⊗ is the Kronecker product. C g and R g are the polarization covariance matrix and interference covariance matrix of ground scattering; C v and R v are the polarization covariance matrix and interference covariance matrix of canopy scattering [29]. In data processing, the MBMP covariance matrix is estimated with the samples in L multilooks:Ŵ Then, considering the singular value decomposition (SVD) ofŴ after rearrangement [29], the MPMB covariance matrix can be decomposed aŝ where P represents the number of scattering mechanisms. λ p is the corresponding singular value. U p and V p are obtained by reshaping the singular value vectors u p and v p into matrices 3 × 3 and N × N elements, respectively. In this particular work, P = 2 because only the signal contributions from the ground and canopy are considered. Next, sort λ p in descending order and define R p = conj(V p ); the interference covariance matrices of ground scattering and canopy scattering can be solved through a combination of linear equations with P = 2: where a and b are two real numbers. Finally, a and b can be optimized using the model-free separation under the constraint that both R g and R c are semipositive definite

SKP-IMLE Pol-TomoSAR Estimator
After the above SKP decomposition, we reshape the observation signals according to (3). Each of them contains only one kind of specific scattering mechanism, i.e., where G g represents the observation data set containing only the surface scattering mechanism, and G c represents the observation data set containing only the canopy scattering mechanism. γ g and γ c are the corresponding reflectivity function; e g and e c are the corresponding noise. It follows after (16) that the covariance matrices of different scattering mechanisms can be obtained: where K g and K c are the backscattering power matrices of ground scattering and canopy scattering. δ 2 g and δ 2 c are the noise variances of ground scattering and canopy scattering, respectively. In this particular work, the decomposition error is ignored.R g andR c are approximated by the two covariance matrices in (14), respectively.
Substituting the covariance matrix of different scattering mechanisms in (14) into the IMLE algorithm and the TomoSAR inversion result based on the SKP-IMLE method can be obtained.
If the tomograms of different scattering mechanisms are obtained using the SKP-IMLE method, the corresponding phase centers can be extracted from them, and then the inversion results of forest height and underlying topography can be obtained.
More specifically, the way to estimate a digital terrain model (DTM) is as follows: For the estimation of the canopy height model (CHM), an additional operation is required because the top position of the tree is higher than the position of the canopy scattering center: Consider the method in [11]. The canopy phase center height h c is calibrated by performing a fitting analysis of the top height h top and the phase center height h c in some samples.
where m and n are the coefficients of the equation. Finally, subtract h g from h top to get the estimated value of the forest height (CHM). In order to explain the proposed TomoSAR method more clearly, the specific process of SKP-IMLE spectral estimator is shown in Table 1. Table 1. Details of the sum of Kronecker products (SKP) decomposition and iterative maximum likelihood estimation (SKP-IMLE) estimator. TomoSAR: synthetic aperture radar tomography.

Study Area and Dataset
In order to verify the correctness and feasibility of the SKP-IMLE Pol-TomoSAR method, the P-band full polarimetric focused SAR images of the TropiSAR2009 project were used for the real data experiments. The study area was located in Paracou, French Guiana, and its geographic location is shown in Figure 1. This test area was a tropical forest containing many tree species. The forest height in the test area was roughly between 20 m and 45 m, and the altitude was between 5 m and 50 m. The test area has six passes of fully polarimetric SAR images that were acquired by the ONERA SEFHI airborne system. The baseline length between adjacent orbits was about 15 m, and the incident angle variety was between 20 degrees and 60 degrees. The system parameters and baseline parameters of the SAR data in the test area are shown in Table 2. In addition, the French Agricultural Research Center for International Development (CIRAD) enabled the light detection and ranging (LiDAR) data to cover a small part of the test area, which was used to verify the correctness and accuracy of the inversion results. In order to evaluate the accuracy of the forest height and underlying topography retrieved by the proposed method, two test lines and two regions of interest (ROI) with different mean forest height were selected for inversion. The mean tree height of ROI1 was 28.9 m and the mean tree height of ROI2 was 25.6 m. The positions of the two ROI in the SAR images and the corresponding averaged LiDAR DTM and CHM are shown in Figures 2 and 3.

Tomograms of the Selected Azimuth Profiles
In order to verify the high resolution of the IMLE method in single-polarimetric TomoSAR and its ability to focus the canopy phase centers and the ground phase centers, two azimuth profiles (the red dashed lines in Figure 3) were selected for estimating the tomograms of HH channel and HV channel, respectively.
In the inversion process, the window size used to estimate the covariance matrix of multibaseline SAR data was 39 × 39 pixels, which roughly corresponds to the 50 m × 50 m in the ground range and azimuth plane. In order to keep the balance between inversion accuracy and calculation efficiency, the maximum number of iterations was set to 10. Figure 4 shows the normalized tomograms of the two test lines at different locations estimated using the IMLE method. All the tomograms were clear and the phase centers of different scattering mechanisms were clearly distinguished. It can be seen from the tomograms that the scattering power of the HH polarimetric channel was mainly concentrated on the ground (as shown in Figure 4a,c), while the scattering power of the HV polarimetric channel was mainly concentrated on the canopy (as shown in Figure 4b,d). In addition, the tomograms of all polarimetric channels were in agreement with the LiDAR data. The inversion results of two test lines at different positions also confirmed that the IMLE method had good reconstruction performance at both near slant-range and far slantrange locations. The continuous and highly focused phase centers of different scattering mechanisms were obtained from both polarimetric channels (see Figure 4). However, due to the influence of forest density, electromagnetic wave power, forest height, and other factors, there was no guarantee that vertical profile could generate two clear and separate peaks for all resolution cells.
In order to obtain more accurate inversion results of forest height and underlying topography, as well as to make the inversion strategy more stable, the SKP-IMLE Pol-TomoSAR estimator proposed in this paper could be used to estimate the tomograms of different scattering mechanisms for the two test lines. Figure 5 shows that the SKP-IMLE Pol-TomoSAR method obtained more focused tomograms than the IMLE method for a single polarimetric channel. The tomograms of canopy scattering and ground scattering at two test lines were continuous, which was consistent with the LiDAR data. The phase centers of different scattering mechanisms were accurately located from the corresponding tomograms. Most importantly, the SKP-IMLE method ensured that we extracted the corresponding phase centers from the tomograms in all resolution cells.

Forest Height and Underlying Topography Inversion
Based on the high-quality tomograms obtained using the SKP-IMLE method, the phase centers of different scattering mechanisms were extracted. Then, the underlying topography could be inverted based on the height values of the ground phase centers, and the forest height could be inverted based on the height values of the canopy phase centers. It should be noted that the locations of treetop were above the locations of the canopy phase centers (as shown in Figure 5b,d). The extracted canopy phase center positions still need to be corrected. For this purpose, 18 samples were selected in the test area, and the average height of the LiDAR-based digital surface model (DSM) in each sample was used to calibrate the phase center height according to (20). After that, the DSM was inverted via the SKP-IMLE method was obtained according to the calibration criterion (as shown in Figure 6). The terrain height was subtracted from the DSM to obtain the estimated value of the forest height.

Comparison Experiments
In order to compare the reconstruction performance of the proposed method and the traditional Pol-TomoSAR methods (such as Pol-beamforming and Pol-Capon), the retrieved tomograms of test line1 were selected for analysis, and the forest height and underlying topography were extracted from the tomograms. The specific processes of the polarimetric Capon and beamforming methods were previously introduced in [24,25]. It should be noted that the window size used to estimate the covariance matrices in the beamforming method and Capon method was the same as that of the SKP-IMLE method.
Taking the tomograms of test line1 as an example (as shown in Figure 9a,c), SKP-IMLE obtained continuous and accurate phase centers. At the same time, the forest height and underlying topography retrieved from the tomograms were in agreement with the LiDAR measurements.
Although the Pol-beamforming method and the Pol-Capon method have the ability to distinguish scatterers with different scattering mechanisms, they still cannot achieve accurate positioning of the phase centers. In some resolution cells, only the phase centers of one scattering mechanism can be detected, and the inverted phase centers are not continuous (see the blue elliptical areas in Figure 10a,c). Therefore, the underlying topography and forest height information extracted by these two methods were not continuous (see the blue elliptical areas in Figure 10b,d).   Figures 11 and 12 show the relationships between the inversion results and the LiDAR measurements. In addition, the root-mean-square error (RMSE) between the results and the LiDAR measurements is shown in Table 3. Figure 11 shows the results of three TomoSAR methods in ROI1. At the same time, the correlation coefficients between the inversion results and the LiDAR measurements were also calculated. The experimental results show that the inversion results of SKP-IMLE were the most relevant to LiDAR measurements (as shown in Figure 11e,f), and had the highest inversion accuracy (as shown in Table 3).
The mean tree height in ROI2 is lower than that of RO11, and the canopy phase centers may be closer to the ground phase centers. Therefore, the inversion of forest height and underlying topography in ROI2 requires higher vertical resolution. The experimental results show that the SKP-IMLE method maintained a good performance in this area (as shown in Figure 12e,f). However, the reconstruction performances of the traditional methods were significantly reduced. In particular, the inversion result of forest height became very unreliable.

Discussion
In the experiment using TropiSAR2009 SAR data set, two different ROI and two different profiles were chosen to verify and analyze the proposed SKP-IMLE TomoSAR method. The experimental results showed that the new method obtained continuous tomograms of different scattering mechanisms in forest areas (as shown in Figure 5). Based on high-quality tomograms, accurate forest height and underlying topography could be retrieved (as shown in Table 3). The performance of the proposed method was also compared with other methods. We found that the main factor affecting the inversion accuracy is the resolution of the TomoSAR method. Due to the superior performance of SKP-IMLE, it obtained higher inversion accuracy in both ROI. In ROI1, the RMSE of the estimations from the SKP-IMLE method was 1.489 m for the underlying topography and 1.997 m for the forest height, which were both much less than the results of traditional spectral analysis methods. In ROI2, the inversion performances of the traditional methods were seriously degraded, but the SKP-IMLE method still obtained high-quality inversion results (as shown in Figure 12). In ROI2, the RMSE of the estimations from the SKP-IMLE method was 1.786 m for the underlying topography and 1.765 m for the forest height. Due to the limitation of vertical resolution, the traditional TomoSAR methods had low accuracy in the inversion of forest height and underlying topography. The underlying topography was obviously overestimated, while the forest height was obviously underestimated (as shown in Figures 11 and 12). It was usually difficult to identify two separate phase centers in all pixels, resulting in large errors in the positioning of the phase centers.
Although the proposed method has many advantages, there are still some limitations in actual forest monitoring. The stability of the algorithm under different baseline configurations needs to be discussed because there are several matrix-vector operations, including matrix inversions. In addition, the algorithm needs to be optimized to improve its computational efficiency and improve its applicability in large-scale forest height inversion.

Conclusions
In this paper, we proposed a Pol-TomoSAR method based on SKP decomposition and maximum likelihood estimation. The new method was used in the inversion of forest height and underlying topography over a tropical forest, and reliable results were obtained. The performance of the proposed method was compared with two traditional TomoSAR methods. Through real data and comparative experiments, the following conclusions can be drawn: (1) The SKP-IMLE method has a high vertical resolution, so it can clearly distinguish the scatterers with different scattering mechanisms and achieve high focus of canopy scattering and ground scattering, which is very beneficial for forest height and underlying topography extraction. However, the traditional TomoSAR methods (such as beamforming and Capon) have limited resolution, and there is no guarantee that the vertical profile generates two clear and separate peaks for all resolution cells. These shortcomings make the extraction of forest height and underlying topography using traditional spectral analysis methods prone to errors. (2) The proposed SKP-IMLE method can obtain good forest height and underlying topography inversion results in tropical forests. With respect to LiDAR measurements, the SKP-IMLE method achieved high inversion accuracy in the two regions of interest.
In conclusion, the SKP-IMLE method has a high vertical resolution and can obtain reliable forest height and underlying topography inversion results even in dense tropical forest areas. Finally, further work will focus on using SKP-IMLE in SAR tomography with a small number of measurements. In addition, the proposed method will be used to estimate the height of boreal forests to test its applicability in different types of forests.