An improved adaptive decomposition method for forest parameters estimation using polarimetric SAR interferometry image

ABSTRACT Forest parameters estimation using polarimetric synthetic aperture radar interferometry (PolInSAR) images is one of the greatest interests in remote sensing applications. Applying the model-based decomposition concept to PolInSAR data opened a new way for forest parameters estimation. However, the method tends to underestimate the forest height due to reflection symmetry assumption and required the numerical solution of nonlinear equation system. In order to overcome these limitations, an improved adaptive decomposition technique with PolInSAR data is proposed. In this approach, an accurate topographical phase and asymmetry volume scattering model are applied to the model-based decomposition technique for polarimetric SAR interferometry image. The accurate topographical phase is first estimated and then used as the initial input parameter to our numerical method. This approach is not only avoiding large error generated by the constant topographical phase in fluctuating forest areas but also improve the accuracy of forest height estimation and the magnitude associated with each mechanism. The performance of proposed method is demonstrated with simulated data from PolSARproSim software and SIR-C/X-SAR spaceborne PolInSAR images over the Tien-Shan areas, China. Experimental results indicate that forest parameters could be effectively extracted by proposed method.


Introduction
The estimation of forest parameters and vertical forest structures plays in an important role in the forest management at local scales and in climate modeling at regional and global scales. Also, the results of estimation contributes scientific data to research of the carbon cycle at region and local scales. In recent years, a new active microwave remote sensing technique based on polarimetric synthetic aperture radar interferometry (PolInSAR) systems has proved to be an effective tool for retrieving forest heights because of its sensitivity to the vertical structure and physical characteristic of the scattering media. In two recent decades, a great number of approaches of forest height estimation using single baseline PolInSAR data have been proposed. These approaches can be divided into two main categories: inversion processing and model-based decomposition. The first group is based on forest heights inversion approaches as introduce by Cloude (Cloude & Papathanassiou, 2003), Yamada (Yamada, Yamaguchi, Kim, Rodriguez, & Boener, 2001a, 2001b, Yamada, Yamazaki, & Yamaguchi, 2006, and Garestier (Garestier & Toan, 2010). The inversion approaches are based on the sensitivity of the radar interferometric measurement to the vertical distribution of scattering elements, combined to the sensitivity of radar polarimetry to shape and orientation of these elements. To invert the PolInSAR measurements into forest height, a forward coherent model has been developed, assuming canopy layer as a homogeneous volume constituted by randomly oriented particles in the whole height interval. These methods are quite simple and most widely used. However, these methods tend to underestimate attenuation of electromagnetic in the ground medium and accuracy of these method become inappropriate for dense forest region due to overestimated of the volume scattering contribution.
The second major group is based on decomposition techniques for PolInSAR. The ideal of applying the model-based decomposition to PolInSAR data was first proposed by Bermand (Ballester-Bermand & Lopez-Sanchez, 2010). It was assumed that the cross-correlation matrix obtained from PolInSAR observations may be expressed as a sum of three scattering mechanism matrices, each representing the contribution of single-bounce, double-bounce and volume scattering. This algorithm enables the retrieval of not only the power contribution from three scattering mechanisms but also the determining their locations in the vertical dimension. However, this method used an assumption of flat topography (i.e. a constant topographical phase) which cannot be suitable for forest areas with relief terrains. Because the decomposition results such as forest heights directly depend on topographical phases (Wilsen, Sarabandi, & Lin, 1998). Therefore, decomposition of PolInSAR measurement data of fluctuant forest areas using this assumption can produces significant errors. In addition, this approach does not utilize full polarimetric interferometry information, so the numerical solution of a nonlinear equation system is required for successfully employing the model-based decomposition algorithm. Another limitation of this algorithm is that an assumption of scattering reflection symmetry held during the observation for all volume scattering terms.
We has recently proposed an adaptive model-based decomposition for PolInSAR data (Minh, Zou, & Yan, 2012). In this method, we suggest that the reference volume scattering covariance can be used to determine the best fit parameters to express general volume scattering covariance matrix. However, there exit an underestimated problem in these methods. One parameter in the surface or double bounce scattering model is set to zero, thus leading to the instability of the decomposition. The instability of these method will cause false forest parameters estimation.
For these reasons, we propose an improved adaptive model-based decomposition technique for forest height estimation using L-band PolInSAR data, where no reflection symmetry required. To achieve this ultimate goal, we use a generalized volume scattering model that is introduced by Arri (2011;Arri, VanZyl, & Kim, 2010). The model does not require any geophysical media symmetry assumption but in order to estimate parameters of the model somewhat complicated computation is required. The cross-covariance matrix obtained from PolInSAR observables is decomposed into the three scattering elements proposed by Freeman and Durden for polarimetric SAR (PolSAR) data. The proposed method is executed following four steps. Firstly, the actual topographical phase is estimated by total line-fit square technique, which is used as initial parameter of our numerical adaptive decomposition method. Secondly, the parameters of three scattering mechanisms are estimated by using adaptive decomposition technique. In this step, the topographical phase is used to solve the nonlinear equation system for target decomposition by using Newton-Raphson method. Thirdly, the forest height is estimated by phase differencing between the canopy phase and the actual topographical phase. Finally, the forest height is compensated by the coherence amplitude approach. The proposed method not only enable the enhancement the accuracy of forest parameters but also of the magnitude associated with each scattering mechanism. Another advantage of the proposed method is that it allows a more robust implementation and an unambiguous estimation of the ground topography as well as canopy phase. Additionally, the proposed method avoids large error generated by the constant topographical phase in fluctuating forest areas and fully utilized the polarimetric interferometry information contained in measured PolInSAR data.
The organization of this paper is as follow. Section 2 introduces the method use to estimate the topographical phase and the improved adaptive decomposition technique for estimation of forest parameters. The experimental results of the proposed method with simulated data and space-borne data are presented and discussed in section 3. Finally, the conclusion and future work are drawn in section 4.

The improved adaptive model-based decomposition
Topographical phase estimation using total least square line fit method A fully polarimetric interferometry system measures for each resolution element in the scene from two slightly different look angles, two scattering matrices S 1 ½ and S 2 ½ . In the case of backscattering in a reciprocal medium, the 3D lexicographic scattering vectorsk L 1 andk L 2 are given by (Cloude & Papathanassiou, 1998) Where i ¼ 1; 2 denote the master and slave images, respectively. The superscript Á ð Þ T represents the matrix transpose and S iHH ; S iHV ; S iVV are elements of the scattering matrices S i ½ , which are also referred to as the complex scattering coefficients.
The complex information measured by the PolInSAR system can be presented in form of three 3 × 3 complex matrices C 11 ½ , C 12 ½ , and C 22 ½ formed using outer products ofk L 1 vàk L 2 as Where : h i denotes the ensemble average in the data processing, and Á ð Þ ÃT denotes the matrix conjugation transpose. Two matrices C 11 ½ and C 22 ½ are the conventional Hermitian covariance matrices that describe the polarimetric properties for each image separately, while C 12 ½ is a non-Hermitian complex matrix that contains polarimetric and interferometric information. The baseline-copolar coherence function is expressed follow (Cloude & Papathanassiou, 1998): Whereω i ; i ¼ 1; 2 f g are the unitary complex vector that defines the polarization of the master and slave images. According to Equation 3, for any polarization vector, the coherence can be calculated from the covariance matrix as a ratio of their quadratic forms. The range of this function is called the coherence region (Cloude, 2009). The ground topographic phase is estimated by performing a linear fit along axis of the coherence region.
In the total line fit square (TLS) method, Cloude used multiple polarization channels and a least square line fitting the multiple complex data points. This is a more robust version of line fit techniques (Cloude, 2009) because it allows to avoid problems of any selected pair of points too close, and thus minimizes the surface topography estimation errors. For ensuring linearity between the real and imaginary parts of coherence, the linear coherence loci assumption is used, as shown in Figure 1.
Errors of two estimated coefficients M and C of line fit can be reduced by implementing a total least squares solution that accounts for errors in both x and y. Geometrically, the TLS approach amounts to using a different measure of distance: the perpendicular distanceR i , at right Figure 1 and related to Δy as shown in Equation 4.
After that, we can be obtained values for the estimate of M and C by minimizing the function in Equation 4. This provides us with a method for fitting a line to an arbitrary number of polarization channels. Finally, we use the estimates of M and C to find the two-unit circle intersection points as show in Figure 2. These two points can be found explicitly in terms of M and C which is given by: The TLS approach employs several typical polarizations, such as HH, VV, HV, HH-VV, HH +VV. . . to do the total least squares line fit for the extraction of the line segment mentioned in the forest model. This approach is often combined with the coherence separation optimization method to decide which phase point to use as the topography estimate. In the proposed method, the optimal coherence coefficients γ i opt i ¼ 1; 2; 3 f g of the PolInSAR system are instead of the complex coherence coefficients γ VV ; γ HHþVV ; γ HHÀVV for the channels VV, HH + VV, HH -VV to improve accuracy when determining the best-fit straight line corresponding to the complex coherence coefficients. These optimal coherence coefficients correspond to the square root of eigenvalues of the matrix K ¼ T À1 22 Ω Ã 12 T À1 11 Ω 12 (Zhu, Zhang, & Li, 2016). Therefore, we use the complex coherence coefficients γ i opt i ¼ 1; 2; 3 f g and γ HV ; γ HH to find out the topographic ground phase.

Improved adaptive model-based decomposition of PolInSAR data
As proposed by Ballester-Bermand and Lopez-Sanchez (Ballester-Bermand & Lopez-Sanchez, 2010), the cross-correlation matrix C 12 ½ , which includes both polarimetric and interferometric information, can be decomposed into the three matrices C S ½ , C D ½ and C V ½ , which correspond to the single bounce, double bounce, and volume scattering mechanisms, respectively.
Subsequently, the three scattering mechanisms to polarimetric interferometry observables are analyzed.  Note that, as already demonstrated by Bermand (Ballester-Bermand & Lopez-Sanchez, 2010), a null correlation between copolar and cross-polar channels will be assumed.
The single bounce scattering model is presented by the first-order Bragg surface scatter, plate. Sphere scattering modeling is slightly rough surface scattering in which the cross-polarized component is negligible. The scattering coefficient amplitude is unchanged for the both images, except the difference in the phase term. This phase term has two contributions: the difference due to the complex scattering coefficient on the case of using different polarizations for master and slave images φ HV ¼ φ H À φ V , and the interferometric phase related to the position in the vertical coordinate ϕ S . Hence, the covariance matrix for the single bounce contribution is expressed as (Guo, Li, Zhang, Yin, & Hong, 2015): In this case, the surface scattering is modelled under the Bragg condition (Freeman & Durden, 1998). Additionally, the Bragg model states that VV j j> HH j j, so these scattering parameters for master and slave images depend only on the surface complex permittivity and the incidence angle. There is the same surface complex permittivity for the master and slave images, and only a slight difference between incidence angles. Therefore, the scattering coefficient amplitude is unchanged for both images, except the difference in the phase term. Thus, it can be assumed that S 1HH j j¼ S 2HH j j¼ S HH j j, where ϕ S denotes the interferometric phase.
The double bounce scattering component is modeled by scattering from the interaction between the ground and tree trunk, where the reflector surface can be made of different dielectric materials. Hence, the covariance matrix of the double-bounce scattering contribution is given by (Guo et al., 2015): Where F D ¼ R GV R TV j j 2 e jϕ D and α ¼ R GH R TH R GV R TV e jðφ V Àφ HÞ . The phase ϕ D denotes the interferometric phase for the double bounce scattering component, and φ V À φ H is the phase difference between the different polarization channels. The coefficients R GH and R GV are the horizontal and vertical Fresnel reflection coefficients of the ground surface, respectively. Similarly, the vertical trunk surface has reflection coefficients R TH and R TV for horizontal and vertical polarization, respectively. These coefficients assumed to be equal for both ends of the baseline.
The volume scattering is direct diffuse scattering from the canopy layer of forest model. The scattering from the forest canopy layer can be theoretically characterized by a cloud of randomly oriented infinitely thin cylinder, and a uniform probability function for orientation angle (Freeman & Durden, 1998). However, for forest areas where the vertical structure seems to be rather dominant, the scattering from tree-trunks and branches display a nonuniform angle distribution. Consequently, the volume scattering contribution is assumed as the nth power cosine-square distribution of orientation with probability density function as introduced by Arri (Arri et al., 2010). This function can be described by two parameters: the particle mean orientation θ, and orientation randomness degree ν. The former θ 2 0 Ä π=2 ½ , and the latter changes in a range of between 0 and 0.91. In order to improve the general for volume scattering contribution, we add the particle scattering anisotropy into the volume scattering contribution proposed by Arri (Arri et al. 2011). The particle scattering anisotropy is a characteristic of the effective shape of an average particle and depends on the particle and background permittivities. If the particle permittivity is similar to the background permittivity, ε r % 1, the particle scattering anisotropy tends toward zero independently of the real shape of the particle and the scattering effectively vanishes. In the other case ε r ) 1 ð Þ, assuming simple ellipsoid particles, one can make the following predictions about the effective particle shapes: as δ j j ! 0, the average effective particle shape approaches an isotropic sphere/surface, whereas for δ j j ! 1 the effective shape of the scattering particle in the polarization plane tends toward a dipole (Neumann, Famil, & Reigber, 2010). Therefore, C v is a generalized volume scattering matrix, which depends on the mean orientation angle, degree of randomness, and the particle scattering anisotropy δ. Then, the volume scattering covariance matrix is given by Where the coefficient p v ð Þ and q v ð Þ are characteristic by sixth-order polynomials as introduced by Arri (Arri et al., 2010). The basic coherence matrices C a ; C b , and C c are expressed as: The particle scattering anisotropy magnitude δ j j is directly related Cloude'sαangle (Cloude & Pottier, 1997) hv iÞ. Based upon the asymmetry volume scattering model in Equation 10, we shall develop the approach that can improve the adaptive model-based decomposition of PolInSAR data (Minh et al., 2012). By replacing the volume scattering term in Equation 7 by Equation 10, we can obtain the new adaptive model-based decomposition for PolInSAR image as follows: We first implement finding the volume scattering covariance matrix with each pair of value θ; σ À Á in their entire range. And then, we use the following equation to find the coefficient F V : Bermand (Ballester-Bermand & Lopez-Sanchez, 2010) analytically derived eigenvalue of C 0 remainder ½ to find the maximum F v that can used in Equation 12. However, an analytical tractable solution of their method also depends on the reflection symmetry assumption. In this case of no reflection symmetry, no computationally efficient algorithm has been derived yet. In this paper, no reflection symmetry is assumed so that the full 3 × 3 covariance matrix is considered. Unfortunately, the analytical way to find the maximum F V is no longer straightforward as the characteristic equation of remainder matrix now is a general cubic polynomial. Instead, we calculate the coefficient F V by marking the following assumptions about three component scattering problem: (1) C S ½ and C D ½ are two unknown rank-1 matrices, (2) C V θ; σ À Á Â Ã is a known positive-definite Hermitian matrix at specific randomness and mean orientation angle, and (3) the best fit F V under the condition with the C remainder ½ matrix becomes zero. Consequently, Where det Á ð Þ denotes the matrix determinant. We can show that Equation 14 is a general cubic equation about F V whose three roots can be easily obtained as in Appendix. In fact, the phase center separation of volume scattering contribution can lie anywhere between the halfway and to the top of canopy layer. That mean as the phase center separation increases due to changes in structure function so, at the same time, the effective volume depth decreases (as the structure function becomes more localized near the top of the layer), and hence the level of volume decorrelation will decrease (Cloude, 2006(Cloude, , 2009. Therefore, we can show that the maximum F V equivalent to best fit under condition that the argument associated with the argument term of the complex coefficient F V is maximum.
. When the parameters F v is selected, we can remove the volume scattering component from the original cross covariance matrix. The remainder matrix is presented as: As can be seen, in matrix C 0 remainder ½ there appear four complex unknowns and four complex observables, since C remainderð1;3Þ ÞC remainderð3;1Þ . Equation 16 lead to a determined nonlinear equation system. Therefore, to determine the rest of unknown parameters α; β,F D ,F S , and C remainder ½ simultaneously, a Newton-Raphson algorithm is implemented. Assuming that the phase centers of double-bounce scattering mechanism is located at ground level [5]; therefore, the topographical phase is regarded as the phase of F D . Hence, the double-bounce scattering phase center is used as the initial input parameter to the Newton-Raphson algorithm, and then, it is used to solve the nonlinear equation system for target model-based decomposition techniques. Thus, for each pair of randomness σand mean orientation angle θ 0 , the parameter set F i ; α; β f g i ¼ S; D; V ð Þin this step is determined from condition minimum of Frobenius norm of matrix C remainder ½ ¼ C 12 ½ À P i¼S;D;V F i C i ½ . We show that the optimum parameter set F i ; α; β; θ 0 ; σ; δ Þis equivalent to the best fit under the condition that the Frobenius norm of matrix C remainder ½ becomes zero, where the estimated parameters are perfectly matched to the observations. Finally, we repeat above steps for each pixel in the image. The algorithm is summarized in Figure 3.

Forest height estimation
One of the simplest approaches to forest height estimation is to use the phase difference between interferogram as a direct estimate of height. Based on the obtained optimization parameters from improved adaptive decomposition, the forest height can be extracted by phase differencing between the canopy phase and ground phase, as in Equation 17: Where ϕ V and ϕ 0 are canopy phase and surface phase, respectively. The parameter θis the mean angle of incidence, R is the distance between radar and an observed point, δ is the baseline tilt angle, B is the baseline and λ is the wavelength. In order to improve the accuracy of forest height estimation, we first use forest height, which is estimated by proposed method, exactly as proposed in Equation 16. In 2010, Chen Hong demonstrated that the interferometry phase corresponding to canopy is increase with real forest height (Hong, Hong, & Chao, 2010). The difference between the actual forest and canopy height is called penetrated depth. The penetrated depth depends on the incident angle, forest species, and forest shape. We show that the penetrated depth doesn't change with real forest height. Hence, the scattering center of canopy is always lower than the forest height. This is why we cannot use canopy scattering center to represent the top of the canopy. Hence the true forest height is always underestimated. To progress, one key idea is that this error can be at least partly compensated by employing a coherence amplitude correction term, as introduced by Cloude (Cloude, 2009). Finally, by combining these two terms with a scaling parameter η, we then obtain an approximate algorithm that can compensate variation in structure, as shown in Equation 17 (Cloude, 2009): Whereγ v denotes the complex coherence for the volume alone. In Equation 18, the first term represents the phase coherence while the second term is the coherence amplitude correction. This expression has the right kind of behavior in two important special cases. If the medium has a uniform structure function the first term will give half the height but the second will then also obtain half the true height (if we set η ¼ 0:5). At the other extreme, if the structure function in the volume channel is localized near top of the layer, then phase height will give the true height, and second term will approach zero that reason the weight set as η ¼ 0. To reduce the error from change of extinction coefficient and the vertical structure, we select η ¼ 0:4, as reported by Cloude (Cloude, 2009). Theoretically, the forest height estimation ranges from 0 to 2π=k z when using the proposed decomposition method from L-band PolInSAR data. However, in practice, the forest height estimation accuracy strongly depends on some factors, such as tree density, species, age-related variation, wavenumber and wave extinction in the volume scattering layer. When the extinction decreases, wave interact with a thicker layer of the volume, resulting in a more important volume decorrelation due to an increase of scatterer height diversity, and in a diminution of the phase center height in the volume, until half of its height. Therefore, for proposed decomposition method, there is an extension of the canopy layer from the crown to the ground when the forest height is largely insufficient. Then, the bottom phases of all the layers overlap at the ground, which is unreliable for the proposed decomposition method of forest height estimation. Therefore, the minimum forest height estimation using the proposed decomposition method is approximately a half of the forest height estimated from the amplitude coherence. Only forests higher than 10 m have been investigated to avoid the dominant ground contribution and to ensure non-overlap of bottom phases of the three layers.
This proposed technique produces accurate results because the adaptive decomposition has been applied to every pixel. However, to achieve this ultimate goal, this technique requires slightly complicated computation. This algorithm enables the power contributions from the three scattering mechanisms to be estimated, as well as determining their locations in the vertical dimension. In addition, the generalized volume component can be characterized by three parameters: a mean orientation angle, a degree of randomness and particle scattering anisotropy, so the distribution of vegetation regions is also determined using this method.

Experimental results
In this section, the effective evaluation of the proposed approach is primarily addressed in terms of the retrieved forest height estimation and ground phase. For such a purpose, the effectiveness of the proposed method is first evaluated with the simulated data, which acquired from PolSARProSim software by Mark L. William (William. 2006), Based on this, the method is adjusted to improve its accuracy. After providing reliable results with simulation data, the proposed method is applied to spaceborne data acquired by SIR-C/X-SAR system from National Aeronautics and Space Administration (NASA).

Simulated PolInSAR data
The proposed algorithm has been first tested with a simulated RVoG scenario, named as PINE in the PolSARProSim software at 1.3GHz and at incidence angle of 30 degrees considering difference soil conditions and averaging window sizes. The interferometer is operated at 10 m horizontal and 1 m vertical baseline. The stand height 18 m, and it is located on a 0.1% ground azimuth and 0.2% ground range slope. The forest stand occupies a 0.72,854 ha area and stand density is 1000 stem/Ha. Azimuth and slant range resolution are 1 m and 1.5 m, respectively.
The correlation between two radar signals is the most important parameter of PolInSAR system, which appraises quality of topical interferometry phase and reflect the volume scatter information in PolInSAR. So, it is important for successful implementation of PolInSAR parameter inversion techniques, estimation of forest height, classification of surfaces features and quality of differential interferometry processing. The correlation coefficient between two radar signals is defined by Cloude and Pottier (Cloude & Pottier, 1997): Where s 1 and s 2 are the two radar signals return from PolInSAR system. In practice, for PolInSAR data, the correlation coefficient ρ c can be evaluated by averaging over a small box of pixels in the interferogram. The value of correlation coefficient run between zero and one. The correlation coefficient equal to one indicates that two signals received by the twice antenna pass are relevant completely, otherwise when the correlation coefficient equal to zero the two signals are irrelevant completely. In addition, in the forest height inversion applications, the effect of decorrelation between two radar signals is important, which lead to overestimated height errors with very small value of correlation coefficient. In our experience with forest height estimation applications, the chosen threshold values for range of correlation coefficient is from 0.45 to 0.95 to neglect the effect of decorrelation between two radar signals on the interferometry phase and the inversion of forest height parameters. The average ρ c s 1 hh ; s 2 hh À Á ρ c s 1 hh ; s 2 hh À Á and ρ c s 1 hv ; s 2 hv À Á of simulated PolInSAR data in this section are 0.8723 and 0.8542, respectively. Hence, the decorrelation between two radar signals affects insignificance the interferometry phase and other parameters. So, in this section, we can neglect the effect of decorrelation between two radar signals on the inversion of forest height parameters. Figure 4(a) shows the simulated scenario that it is the coniferous forest zone, with 143 pixels in range and 131 pixels in azimuth. Figure 4(b) shows a red, green, blue (RGB) coding Pauli image of the forest scenario considered and the red line indicates the transection analyzed in this section. The top of image corresponds to far range, which can be identified due to the shadowing effect at the borders of the forest. The forest scenario considered is the placed above a Bragg surface with slightly sloped terrain. Figure 5 is a plot of the forest height estimation of the proposed approach compared with the threestage inversion (Cloude & Papathanassiou, 2003) and adaptive model-based decomposition method (Minh et al., 2012) in the 134 th row of azimuth transect line (the red line in Figure 4(b). Compared with the three-stage inversion and adaptive modelbased decomposition method, the proposed method provides more accuracy results. The three-stage inversion algorithm proposed by Cloude and Papathanassiou (Cloude & Papathanassiou, 2003) is one of the most successful processes for inversion of forestry parameters using PolInSAR image and most widely used because of its simplicity. In this algorithm, the inversion of forest height can be break into three separate stages. For the first two stages, the ground topography phases are retrieved by using the method of line fit with three separate complex coherence values. The forest height and extinction coefficient are estimated based on an observed volume coherence only in the last stage. In this stage, the authors based on assumption that the HV channel is not containing any ground backscattering contribution and volume decorrelation is γ est;v %γ HV exp Àjϕ 0 ð Þ. Then, they established a lookup table of volume coherenceγ v according to two parameters of forest height h v and the extinction coefficient σ. By minimizing the least square difference between the volume decorrelation and LUT, the forest height can be estimated. Therefore, forest height estimation of three-stage inversion process significantly depends on the estimation of model prediction (LUT). For improvement of forest height estimation accuracy, we proposed the adaptive model-based decomposition method for single frequency, single baseline PolInSAR measurement (Minh et al., 2012). In this method, a volume scattering model was performed under reflection symmetry assumption as reference volume scattering model. And then, the volume scattering covariance matrix was computed so that it approximates to the reference volume covariance matrix. However, the accuracy of adaptive model-based decomposition method becomes inappropriate by reflection symmetry assumption. Based on Figure 5 and Table 1, we can  say that the forest height estimation by using proposed method are more accurate and reliable than its by using three-stage inversion and adaptive modelbased decomposition methods.
The forest height estimation by using the proposed method is shown in Figure 6. In this figure it is shown that almost the peak differential of the forest height is located at 18m approximately. The actual forest heights are quite well retrieve, except some pixels are overestimated but the almost of forest height in these pixels all less than 23 m. The real effective forest height will be higher than these values so we can say that these results are acceptable. Figure 6 shows that the forest height estimation by proposed method is located in a range from 14 m to 24m. Likewise, the proposed method provides relative accuracy with small error, and is more accurate for vertical structure variations.
To further evaluate the performance proposed decomposition method in analyzing the effect of polarization channels to scattering mechanisms contributing, 1200 test samples are randomly taken from PolInSAR data. The amplitude contributions of the three scattering mechanisms to the HH and VV channel are presented in Figure 7. From this figure, we show that the amplitude of volume contribution play of role dominant in the both polarization channels. As illustrated in Figure 7(a), for HH channel, the amplitude response is dominated by volume and double bounce scattering while the amplitude of single bounce scattering mechanism is significant low. It is reasonable that the Fresnel reflections at the trunk and ground surface raises the double bounce scattering contribution in HH polarization channel. Especially, the amplitude of double bounce scattering contribution is the highest around some pixels 150, 205, 430, 595, 800, and 1150. This is because that the double bounce scattering mechanism is mainly caused by interaction between trunk, branch, and ground in the forest area. In these pixels, the branches are mainly oriented in the horizontal direction, which raise the amplitude of double bounce scattering contribution. In the VV polarization channel, the volume scattering is the most dominant expect around some pixels 150, 200, 450 600, 800, and 1170, where single bounce scattering exceeded volume scattering. In addition, the double bounce scattering for VV-pol channel is relative low in terms of the Fresnel definition of smooth surfaces within bare surfaces and areas surrounding the volume scattering. This is consistent with the definition of the Fresnel coefficients for the quad-pol data, where the VV contribution is higher than the HH one (Wilsen et al., 1998).
After that, this paper evaluate the impact of vertical wavenumber on the accuracy of the forest height estimation by using proposed method. About 17 groups of simulated data with vertical wavenumber gradient change were taken from PolSARproSim software by changing the angle of incidence, while the rest of parameters remains unchanged. The range of vertical wavenumber is from 0.0433 to 0.4513, which corresponds to the range of the incidence angle from 60°to 20°. And for each group of simulated data, the forest height are estimated by proposed method. The experimental results are represented in Figure 8. As shown in Figure 8, there was a positive correlation between the vertical wavenumber and the forest height inversion error. When the vertical wavenumber range was between 0.06 and 0.10, the forest height inversion error was the smallest, the error was less than 15%, which was the best vertical wavenumber interval. When the vertical wavenumber increased further, the inversion error increased significantly, which was not suitable for forest height inversion. When the vertical wavenumber was small (k z < 0.06), the data independence was bad, and the decorrelation was serious, so the reliability of inversion results decreased. The unfavorable k z values lead to a reduced inversion performance. Accordingly, too low and/or too high k z values have to be masked out during the inversion (in accordance with the expected forest height range).  Figure 6. Forest height is estimated by proposed method.

Space borne PolInSAR data
Next, we have also tested the performance of the proposed method with spaceborne data. The spaceborne data used consists of two SIR-C single look complex (SLC) image pair of the Tien-Shan test site by the SIR-C/X-SAR system on 7 and 9 October 1994 (data takes 122.20 and 154.20). They consist of quad-pol interferometric data at L band with a 24.569-degree angle of incidence and 413.8 m baseline. After coregistration of PolInSAR images, we select evaluation area with 495 pixels in range and 495 pixels in azimuth. We first evaluate the correlation between two PolInSAR images of evaluation area. The average ρ c s 1 hh ; s 2 hh À Á and ρ c s 1 vv ; s 2 vv À Á of evaluated data are 0.7439 and 0.7381, respectively. Hence, the decorrelation between two radar signals insignificantly affects the performance of forest height estimation. The evaluation region has a mixed forestry, road and agricultural area. Figure 9(a) shows the optical image of the test site from Google Earth with the scale, longitude and latitude. Figure 9(b) is a composite image of the evaluation patch in Pauli basis. Figure 10(a) is a plot of forest height estimation of the proposed approach compared with the adaptive model-based decomposition method of a selected row, that is, 150 th row. The parameters of forest used by two methods are calculated and shown in Table 2. From Figure 10, we show that the forest height estimation with the proposed method ranges from 15 m to 27 m, while the forest height estimation using the adaptive method located in a range from 15 m to 22 m. A contrast exits between the results of proposed method and adaptive method in some pixels (such as 40, 50, 98, and 105). A possible reason is that the volume scattering model of the adaptive method tends to the scattering reflection symmetry model for the observation. Moreover, the adaptive method is based on thin cylinders as the primary scattering representing the vegetation canopy. However, many types of vegetation canopy include larger leaves and thicker branches that are not adequately represented by thin cylinder scattering. Whereas, the volume scattering model of proposed method is based on a scattering matrix that could represent more complex scatters and no reflection symmetry is required. Therefore, the forest height estimation by adaptive method is significant lower than its by using proposed method in these pixels. Figure 10(b) shows the estimated forest height by proposed approach in the evaluated area. This figure shows that most of the peak differential of the forest height is located approximately at 20 m. The forest height estimation at some pixels is overestimated but less than 28 m. However, these values are almost lower than the 2π height ambiguities, which about 32 m, so we can say that the results of proposed method are acceptable. Based on Figure 10 and Table 2, we can say that the height forest estimation and underlying ground topographic phase are more accurate and less error prone than its by using the adaptive model-based decomposition method.  The proposed method provides three addition physical parameter maps: randomness, mean orientation, and particle scattering anisotropy. From Table 2, the results of the proposed method indicate that in average the canopy particles can be described as prolonged ellipsoid ( δ j j ¼ 0:7), the mean orientation angle is high ( θ % 1:6), the degree of randomness is relatively high (σ % 0:52). The randomness map is displayed in Figure 11 (a). In the farmland and road areas, the randomness values are close to the delta distribution σ % 0 ð Þ, while these values approximate to cosine-square distribution σ % 0:57 ð Þ in forest areas. It is reasonable seen that the canopy which consists a cloud of randomly oriented scatterers, is more like to induce nonreflection symmetry scattering. Moreover, the radar signal can penetrate foliage and interacts with branches that have higher randomness when compare with its in the agriculture and road areas (Cloude, 2009). The mean orientation angle map is present in Figure 11(b). Pixels with horizontal orientation θ % π=2 À Á are widely distributed in the forest areas. From Figure 11(b), we show that the horizontal orientation angle in the forest area indicates that branches are mainly oriented in the horizontal direction. Pixels with vertical orientation θ % 0 À Á are widely distributed in the agriculture and   road areas. It is because that the single bounce scattering mechanism is strong in the both these two areas. However, mean orientation angle values are around θ % π=2 at some pixels in the agriculture areas. A possible reason that the double bounce scattering mechanism is dominant in these pixels, that it raises the horizontal contribution due to the Fresnel reflections at the trunk and ground surface. Therefore, the horizontal orientation in these pixels is not related to the physical orientation of scatterers in the canopy layer. After that, the magnitude of particle scattering anisotropy map is represented in Figure 11(c). The particle scattering anisotropy is a measure of the effective shape of the particle, as observed in the polarization plane. Theoretically, the magnitude of critical particle scattering anisotropy located in a range from 0 to infinity. However, to make the particle scattering anisotropy result more meaningful for forest areas, we restrict the range of the particle scattering anisotropy value 0 δ j j 1 ð Þ for good pixels that fits the proposed method (Neumann et al., 2010). From Figure 11(c), we show that the magnitude of particle scattering anisotropy δ j j was ranged between 0.5 and 1 in the forest area, whereas it closes to zero in the agriculture and road areas. A possible reason that single bounce scattering contribution is dominant in the agriculture and road regions, which is derived from contribution of spheres or surface scatterers. Moreover, the single bounce scattering model is presented by the first order Bragg surface scatterers, plate, sphere, and triple bounce scattering modeling slightly rough surface scattering in which the magnitude of particle scattering anisotropy in the both these two areas approaches to zero that are associated with low returns from shadow and slightly rough surface. On other hand, the volume scattering contribution is dominant in the forest area, which is modeled as contribution from a cloud of random oriented cylinder-like scatterers. Therefore, the magnitude of particle scattering anisotropy δ j j in this area approximate to 1. Hence, the particle scattering anisotropy of proposed method is adjusted to fit the PolInSAR data. Based on Figure 11 and Table 2, we can say that the forest parameters estimation by using proposed method are more accurate and reliable than its by using adaptive modelbased decomposition method.
Similarly as section 3.1, the results for the amplitude contribution of the three scattering mechanisms are shown in Figure 12(a,b)). For the VV-pol crosscorrelation channel, the volume scattering is the most dominant source of scattering, with direct scattering contributing a little less than volume scattering, except for pixels around 160, 600, and 1100 according to agricultural area or road area, see Figure 12(b). Around these pixels, the volume scattering contribution significant drops due to the increasing of the single bounce mechanism. This is because that the single bounce scattering is dominant in the both these two areas. The amplitude contribution of the double bounce scattering has all low value for the VV channel. However, the double bounce scattering contribution increases and becomes dominant in some areas (such as at some around pixels 500, 510, and 590), which correspond to bare ground (including road and big gaps between trees). It is reasonable since that the Fresnel reflection at trunk and ground raises the double bounce scattering contribution. For the HH polarization cross correlation channel, the amplitude of the double bounce scattering mechanism is relatively high but it is not dominant due to the relative roughness of terrain, while the amplitude of the volume scattering mechanism is also dominant in the forest areas. In comparison with the VV channel, double bounce scattering increased to a level only below that measured for volume scattering. This is consistent with the definition of Fresnel coefficients, in which the HH contribution exceeds the VV contribution.

Conclusion
A method for enhancing the forest height estimation accuracy based on the improved adaptive decomposition technique for PolInSAR data has been proposed in this paper. The method of total least square line fit has been used to estimate the accurate topographical phase of forested areas. This phase is then used as the initial input parameter of the Newton-Raphson method. This method can improve not only the accuracy of forest height estimation but also the accuracy of the fit of the decomposition relative to the physical scattering model. In addition, in forest areas, the oriented angles of scattering from tree trunk and branches exhibit nth power cosine-square distribution and not uniform distribution. Our decomposition is implemented by adjusting the generalized volume scattering model proposed by Arri to produce the close fit to observed data. As a result, the power contribution and phase centers of three scattering mechanism are estimated. The proposed algorithm has been tested using simulated and spaceborne data, and its performance has been evaluated in comparison with the threestage inversion and adaptive decomposition method. The estimated height of the volume scattering obtained from the accurate topographical phase corresponds more closely to the tree height than from the constant topographical phase. Experimental results indicated that vegetation parameters can be retrieved directly and accurately. In the future, we continue to implement more experiments using various data of different fields for the full performance evaluation of the proposed algorithm.

Disclosure statement
No potential conflict of interest was reported by the authors.