Forest height estimation from mountain forest areas using general model–based decomposition for polarimetric interferometric synthetic aperture radar images

Abstract The estimation of forest parameters over mountain forest areas using polarimetric interferometric synthetic aperture radar (PolInSAR) images is one of the greatest interests in remote sensing applications. For mountain forest areas, scattering mechanisms are strongly affected by the ground topography variations. Most of the previous studies in modeling microwave backscattering signatures of forest area have been carried out over relatively flat areas. Therefore, a new algorithm for the forest height estimation from mountain forest areas using the general model–based decomposition (GMBD) for PolInSAR image is proposed. This algorithm enables the retrieval of not only the forest parameters, but also the magnitude associated with each mechanism. In addition, general double- and single-bounce scattering models are proposed to fit for the cross-polarization and off-diagonal term by separating their independent orientation angle, which remains unachieved in the previous model-based decompositions. The efficiency of the proposed approach is demonstrated with simulated data from PolSARProSim software and ALOS-PALSAR spaceborne PolInSAR datasets over the Kalimantan areas, Indonesia. Experimental results indicate that forest height could be effectively estimated by GMBD.


Introduction
Polarimetric interferometric synthetic aperture radar (PolInSAR) is a type of remote sensing technique that integrates the advantage of both polarimetric SAR and interferometric SAR. Therefore, retrieving parameters related to forest and monitoring the underlying soil surface have become one of the major applications in the field of active microwave remote sensing including SAR system. In the last 15 years, many techniques have been proposed for forest height estimation using single-baseline PolInSAR images, most of which can be broken up into two categories. The first group is based on the random volume over ground model as introduced by Cloude and Papathanassiou, 1 Yamaguchi et al., 2 Garestier and Le Toan. [3][4][5] However, these methods tend to underestimate forest height due to attenuation of electromagnetic (EM) waves in the ground medium, and the accuracy of these methods becomes worse for dense forest areas due to the strong volume scattering component. The second major group is based on model-based decomposition technique for PolInSAR image as introduced by Ballester-Berman and Lopez-Sanchez 6 and Neuman et al. 7 The technique opened a new way for forest height estimation using the model-based decomposition technique. In this technique, they assumed that the volume contribution is a cloud of uniformly randomly distributed thin cylinders. However, natural forest canopy shows preferential orientation of the branches, smaller twigs, and leaves, which limits the applicability of the uniform randomly distributed model. 8 So, this model is suitable for forest areas that assumes scattering reflection symmetry, but cannot be generalized for forest areas. One facet of this article is to provide an approach for canopy parameter estimation that is not constrained by these limits.
On the other hand, most previous studies were carried out over relatively flat areas. For mountain forest areas, scattered signals are strongly affected by the variations of the local incidence angle and the local orientation angle due to the local topographic slope. Results of VanZyl 9 and Lin and Sarabandi 10 indicated that the dominant scattering is changed drastically when the ground surface is tilted relative to the horizontal. In Ref. 11, Lee and Anisworth investigated the effect of ground topography on the coherency matrix and polarimetric target decomposition. However, in this method, the double-and single-bounce scatterings are still modeled with zero cross-polarization power. In order to overcome these limitations, general double-and single-bounce scattering models will be proposed by separating their independent orientation angles.
This article aims to present an accuracy improvement of forest height estimation over mountain forest that accounts for the effect of slope terrain on scattering mechanisms. In this article, we proposed a forest height estimation approach from mountain forest areas using general model-based decomposition (GMBD) technique for PolInSAR image. In this approach, we shall develop a general volume scattering model that can be characterized by two parameters: a degree of randomness and a mean orientation angle. 8 After that, we will form a lookup database including reference volume scattering coherence matrices. The general volume scattering coherence matrix can be determined during the optimization procedure of the model inversion. The unknown parameters of forest can be obtained by using the nonlinear least square optimization method. Finally, the forest height is estimated by phase differencing between canopy phase and underlying ground phase. This algorithm enables the retrieval of not only the forest parameters above the tilted ground plane, but also the magnitude associated with each mechanism. In addition, the negative power can be avoided, and the proposed approach allows a more robust implementation and an unambiguous estimation of ground topography. Another advantage of the proposed algorithm is that it makes use of all the information provided by the interferometric coherence matrix, which remains unachieved in the previous model-based decompositions. Experimental results show that the accuracy of the forest height can be improved by the proposed approach.
The organization of this article is as follows. In Sec. 2, the forest scattering model for a sloping forest area is presented. The retrieval of forest height using GMBD algorithm is delivered in Sec. 3. The experimental results of the parametric inversion with simulated data are presented and discussed in Sec. 4. Summary and conclusions are drawn in Sec. 5.

Polarimetric Interferometric Coherence Matrix in Mountain Forest Areas
A fully PolInSAR system is measured for each resolution cell in the scene from two slightly different look angles and two scattering matrices, [S 1 ] and [S 2 ]. For the case of reciprocal medium and monostatic backscattering, the individual polarimetric datasets may be expressed by means of the Pauli target vector 12 where ½ T represents the vector transposition, S pq (p; q ¼ fh; vg) is the complex scattering coefficient, and i ¼ 1; 2 denotes the measurement at two ends of the baseline. The basic radar observable in PolInSAR is a six-dimension complex matrix of a pixel in each resolution element in the scene, defined as shown in Eq. (2) where h·i denotes the ensemble average in the data processing, and * represents the complex conjugation. ½T 1 and ½T 2 are the conventional Hermitian polarimetric coherence matrices which describe the polarimetric properties for each individual image separately, whereas ½Ω is a non-Hermitian complex matrix which contains polarimetric and interferometric information. For a tilted plane, the horizontal vector is no longer parallel to the surface; so, most of the polarization channels (HH, HV, VH, and VV) are affected by the tilted slope. The amount of slope-induced shift in the local orientation angle can be visualized as the rotation of the verticalhorizontal basis vector about the line-of-sight, so that the horizontal vector is again parallel to the surface. The local orientation angle is geometrically related to topographic slopes and radar look angles, and it is a function of the azimuth slope, the range slope, and the incidence angle of the flat terrain. 11 A schematic diagram depicting the radar image geometry is shown in Fig. 1. The local orientation angle change due to the azimuth slope effect can be expressed as where γ and ω are the local ground range and azimuth slopes, respectively. The angle θ 0 is the incidence angle of the flat terrain. For mountain terrain, the variation of the local incidence angle and the local orientation angle due to the local topography will lead to changes in the scattered signal. 9 The polarimetric interferometric coherence matrix and coherence matrix in sloping terrain can be obtained by rotation of an orientation polarization angle. 11 ½ΩðψÞ ¼ ½RðψÞ½Ω½RðψÞ ÃT ¼ where RðψÞ is the unitary rotation matrix for the coherence matrix, and similar rotation transform was presented as in Ref. 11. Fig. 1 Schematic diagram of the imaging geometry for sloping ground surface.

Single-Bounce Contribution from Sloping Surface
In case of the flat surface, the single-bounce scattering model is presented by the first-order Bragg surface scatter, plate, and sphere, and triple-bounce scattering modeling is slightly rough surface scattering in which the cross-polarized component is negligible. The amplitude of the scattering coefficient does not change for both images, except the difference in the phase term. This phase term will have two contributions: the difference due to the complex scattering coefficient in 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 Pauli target vector of single-bounce contribution for both ends of the baseline is represented as cos α s e jφ 1 sin α s e jφ 2 0 cos α s e jðφ 1 þϕ s Þ sin α s e jðφ 1 þϕ s Þ 0 where the coefficients φ 1 and φ 2 account for wave propagation processing. We assume that these coefficients for master and slave images of the interferometric pair are the same. The real coefficient α s ≤ π∕4 depends on dielectric constant and angle of incidence and is given as 13 The phase ϕ s denotes the interferometric phase for the direct scattering from surface. The coefficients R H and R V are the horizontal and vertical Fresnel reflection coefficients of the ground surface, respectively. In this case, the coherence matrix for the single-bounce contribution is expressed as cos 2 α s cos α s sin α s e jφ s 0 cos α s sin α s e −jφ s sin 2 α s 0 0 0 0 For a rugged terrain surface, the presence of topographic variations induces the local incidence angle θ to be a function of the local terrain slope where the angle θ 0 is the incidence angle of the flat terrain. In addition, the presence of topography variations induces a local coordinate system in accordance with the tilted ground surface. In order to address the orientation of ground scattering, terms of the forest scattering model are superseded by the scattering from the slanted ground plane. In this case, the coherence matrix for the single-bounce contribution is obtained from the rotation of an orientation angle χ. 2 cos 2 α s 2 cos α s sin α s cos 2χe jφ s −2 cos α s sin α s sin 2χe jφ s 2 cos α s sin α s cos 2χe −jφ s 2sin 2 α s cos 2 2χ −sin 2 α s sin 4χ −2 cos α s sin α s sin 2χe −jφ s −sin 2 α s sin 4χ 2 sin 2 α s sin 2 2χ Slope also modifies the apparent alpha parameter α s for the surface backscatter. The range slope component will cause an increase (for slopes away from the radar) or decrease (for slopes toward the radar) in the apparent alpha parameter. In this article, we assume that the range slope is away from the radar.

Double-Bounce Contribution from Sloping Surface
The double-bounce scattering component is modeled by scattering from the interaction between the ground and tree trunk or stem, where the reflector surface can be made of different dielectric materials. Hence, the Pauli target vector of double-bounce contribution for master and slave images can be expressed in the form where α D ≥ π∕4 depends on the angle of incidence and the two dielectric constants of the surface and reflector, which is represented as follows: 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 polarizations, respectively. These coefficients are assumed to be equal for both ends of the baseline. Hence, the coherence matrix for the double-bounce contribution is given by with In the double-bounce scattering component, the backscatter ray path no longer includes the two specular reflection mechanisms for two surfaces over sloping terrain surface. For this reason, the presence of slopes can attenuate the backscatter dihedral return and leave only the direct surface component. The slope tolerance of the double-bounce component depends on many factors such as the height of the vertical scatterer and radar wavelength. Instead of considering this phenomenon on a case-by-case basis, we prefer to model surface as some a priori unknown mixture of singleand double-bounce scattering mechanisms. Hence, the coherence matrix for the double bounce in mountain forest areas is obtained by the rotation of an orientation angle η.

General Volume Scattering Model
The volume scattering is direct diffuse scattering from the canopy layer of forest model. In mountain forest areas, volume scattering component is not much affected by the tilt of the ground surface because trees on a slope grow in alignment with gravity and sunlight. So, the scattering from the canopy layer of forest can be reasonably characterized by a cloud of randomly oriented infinitely thin cylinders. 8,15 In the theoretical modeling of volume scattering, a cloud of randomly oriented dipoles is implemented with a uniform probability function for the orientation angles. However, for forest areas where vertical structure seems to be rather dominant, the scattering from tree trunks and branches displays a nonuniform angle distribution. Therefore, we assume the volume scattering contribution with an n'th power cosine-squared distribution of orientation with probability density function as Ref. 16. The parameter n ranges from zero to infinity. In practice, there is little difference between distributions with values of n larger than about 20 or so. It can be shown that the standard deviation σ can be evaluated in terms of n, as in Ref. 8. The standard deviation σ varies in a range between 0 and 0.91. It changes with the corresponding change in distribution function from delta to uniform distribution function for the vegetation orientation angles. Then, the volume scattering coherence matrix is expressed as whereθ is the mean orientation angle of dipoles. The coefficients pðσÞ and qðσÞ are characterized by sixth-order polynomials, as in Ref. 8. The basic coherence matrices ½T a ; ½T b , and ½T c are expressed as

Forest Height Extraction Based on GMBD
In this section, we propose the algorithm for the estimation of forest height in mountain forest areas using GMBD for PolInSAR image. For PolInSAR data, the polarimetric interferometric coherence matrix after rotation by orientation angle is expressed as follows: ½ΩðψÞ ¼ f s e jϕ s ½T s ðχÞ þ f D e jϕ D ½T D ðηÞ þ f v e jϕ v ½T v ðθ; σÞ þ ½T remainder ; where f S , f D , and f v represent the scattering power coefficient of single bounce, double bounce, and volume scattering, respectively. χ and η are the rotation angles for single-and double-bounce scattering, respectively. T remainder is the residual term which should be minimized after the decomposition. The Frobenius norm of remainder component can be used to determine the best fit parameter to express the measured polarimetric interferometric radar data. Therefore, the optimization criteria are min ∶kT remainder k 2 : In this section, we first find the volume coherence matrix from Eq. (20). In forest area, the backscattering of an EM depends on the shape, size, and orientation of leaves, small branches, and tree trunk. Therefore, we will form the lookup database including the volume scattering mechanism models in Refs. 15 and 17-19 as a reference volume scattering model. These models not only are suitable for geophysical media symmetry, but also satisfy the conditions for geophysical media asymmetry.
In this article, we suggest that each of the reference volume scattering coherence matrices in the lookup database can be used to determine the best-fit parameter to express the general volume scattering matrix. For each of the reference volume scattering coherence matrices in the lookup database, we implemented finding of the volume scattering, so that ½T i v ðθ; σÞ approximates to the reference volume scattering coherence matrix by varying the randomness σ and mean orientation angleθ for the entire range. 8 These parameter sets are equivalent to a best fit under condition that subtraction of general volume scattering coherence matrix and reference volume scattering matrix becomes zero. The optimal volume scattering model can be determined during the optimization procedure of this model. When the set of generalized volume scattering coherence matrices is determined based on the lookup database, we can obtain the set of canopy phase and the coefficient f v as follows: As can be seen, in matrix ½Ω remainder , there appear eight real unknowns ff s ; f D ; ϕ s ; ϕ D ; α; φ; χ; ηg, a matrix ½T i remainder , and nine complex observables, since Ω remainderði;jÞ ≠ Ω remainderðj;iÞ . This formulation leads to a determined nonlinear equation system. Therefore, to determine the rest of the unknown parameters ff s ; f D ; ϕ s ; ϕ D ; α; φ; χ; ηg and ½T i remainder simultaneously, the nonlinear least square optimization method is implemented. 20 In order to improve the quality of this proposed method, the initial values for orientation angle of single bounce κ and double bounce η are all equal to − π 4 (Ref. 10). With each pair value (χ; η), the parameters ff s ; f D ; ϕ s ; ϕ D ; α; φg are determined by employing an singular value decomposition of the Ω remainder . The parameter set ff s ; f D ; ϕ s ; ϕ D ; α; φ; χ; ηg in this step is determined from condition minimum of Frobenius norm of matrix ½T remainder i . We show that the parameter set ff s ; f D ; f v ; ϕ s ; ϕ D ; ϕ v ; α; φ; χ; ηg is equivalent to a best fit under condition that the Frobenius norm of remainder matrix becomes zero, where the estimated parameters are perfectly matched to the observations. Finally, we repeat both the above steps for each pixel in the image. The algorithm is summarized in Fig. 2. When the optimal parameters are obtained from GMBD, we can extract the surface topography phase as in Eq. (27) Based on the obtained optimization parameters from general decomposition approach, the forest height in mountain forest areas can be extracted by using the phase differencing in Eq. (29) where θ is the mean angle of incidence, R is the distance between the radar and an observed point, δ is the baseline tilt angle, B is the baseline, and λ is the wavelength.

Experimental Results and Discussion
In this section, the effective evaluation of the proposed approach is addressed in terms of the retrieved forest height estimation and ground phase. For such a purpose, we have applied the proposed method to a dataset acquired from PolSARProSim software by William 21 as well as spaceborne data acquired by the ALOS-PALSAR system.

Simulated PolInSAR Data
First, the proposed approach has been evaluated by the simulated forest scenario and considering different slope terrains, which are generated with the PolSARProSim software. The simulated data are realized at 1.3 GHz, and the interferometer is operated at 10-m horizontal and 1-m vertical baseline. The stand height 18 m is located on a 18-deg ground azimuth and 6-deg ground range slope. The forest stand occupies a 0.82745 ha area, with a stand density of 360 stem/ha. Azimuth and slant range resolutions are 1.5 and 1.06 m, respectively. Figure 3(a) shows a red, green, blue (RGB) coding Pauli image of the forest scenario considered, and the red line indicates the transection analyzed in this article with 139 × 129 pixels. Figure 3(b) is a plot of the forest height estimation by using GMBD compared with the three-stage inversion and adaptive decomposition methods in the 124th row of azimuth transect line. 1 Table 1 indicates that the proposed method is more accurate and has less error than the threestage inversion and adaptive decomposition methods. The three-stage inversion is the most used coherence model for coherence optimization, and it requires multiple-parameter least square estimation, which is complex and often ill-conditioned. Besides, this method assumed the vertical structure, and minimum ground-to-volume scattering ratio needs to be lower than −10 dB to secure around 10% accuracy. 1 Otherwise, the forest height and extinction estimation by using three-stage inversion process is not reliable, and only the underlying ground topographic phase is reliable. Based on Fig. 3 and Table 1, we can say that the forest height and ground phase estimations by using GMBD are more accurate and reliable than that by using three-stage inversion process and adaptive decomposition.
Changes in the scene parameters can be noticed by means of the proposed method. Table 2 shows the forest parameters estimation with difference slope terrains. The rest of the parameters remain unchanged. The interferometric phase is affected remarkably by both azimuth and ground ranges slopes. So, for the forest height estimation methods related with phase, it will be hard to obtain the right value. 22 In the proposed method, the accuracy of interferometric phase is significantly improved by orientation compensation and by choosing the best-fit parameter for each scattering mechanism using the solving of the nonlinear least square optimization. From Table 2, we show that the local ground range slopes γ from 0 to 30 deg and the local ground azimuth slope ω all range from 0 to 45 deg. Table 2 shows that the forest height decreases when increasing azimuth slope. Especially, the forest height is less accurate at ðγ; ωÞ ¼ ð30 deg; 45 degÞ, but the accuracy of it is relatively high, at 89.96%. Based on Table 2, we can say that the forest height and ground phase estimations by using GMBD are relatively accurate and reliable.

Spaceborne PolInSAR Data
Next, the effective evaluation of the proposed approach is addressed with spaceborne data. The dataset used in this section is acquired from an image pair of the Kalimantan region, Indonesia, by the ALOS-PALSAR system observed on March 12 and April 27, 2007, respectively. The baseline of the two observations is 330 m at the scene center, and the spatial resolution of the test data is 30 × 10 m 2 . They consisted of full polarized data at L-band with 21.5-deg angle of incidence and composed of 12;816 × 1129 pixels. The optical image from Google Earth and the color image of the classical Pauli decomposition are shown in Figs. 4(a) and 4(b), respectively.
The Kalimantan area contains heterogeneous objects such as forest area, agricultural area (violet area), and mountains covered with trees. As analyzed by Papathanassiou and Cloude, 23 the presence of temporal correlation coefficient leads to a decrease of the amplitude of the interferometric coherence, but does not affect the position of the effective phase centers. Furthermore, the amplitude of the interferometric coherence of test site data (after co-registration  Table 2 Forest parameter estimation for difference slope terrains. Parameters True image and filtering procedures) is almost greater than 0.6, and the forest height of the proposed method is estimated by using the difference phase method. Hence, in this section, we neglect the effect of the temporal decorrelation. After co-registration of PolInSAR image, we select two regions of interest from Fig. 4(b) for further comparison including mountain forest area A with 375 × 384 and forest area B with 265 × 481. Patch B contains mainly pure forest, whereas patch A includes mostly mountains covered with trees. Figure 5(a) is a plot of the forest height estimation of the proposed approach compared with adaptive model-based decomposition 24 with orientation angle compensation in patch A. For a fair comparison, all these methods are implemented on the polarimetric interferometric coherence matrix. This figure shows that the forest height of the proposed approach is located in a range from 16 to 24 m, whereas the forest height of the adaptive model-based decomposition method is located in a range from 3 to 23 m. The parameters of forest over two patches are calculated and shown in Table 3. This table indicates that the proposed approach is more accurate and less error prone than the adaptive model-based decomposition approach, even though the adaptive model-based decomposition was compensated by the rotation coherence matrix for the orientation effect. Orientation angle compensation can reduce these cross-polarization powers and better decomposition performance can be achieved. However, the estimated polarization orientation angle is a mixture among all the scattering mechanisms for a coherence matrix. So, this approach cannot always guarantee for cross-polar the single-and double-bounce scattering contributions back to zero. This is a possible reason that the canopy phase and underlying ground phase estimation ambiguities still appear, especially at the mountain forest areas. The results concerning the retrieval of the fraction fill canopy of trees r h are presented in Table 3. The fraction fill canopy of trees is determined as follows: Figure 5(b) shows the estimated forest height by the proposed approach in mountain forest A. This figure shows that most of the peak differentials of the forest height are located approximately at 20 m. The forest height estimation at some pixels is overestimated by less than 25 m. The real effective tree height will be higher than these values, so we can say that the results are  acceptable. Based on Fig. 5 and Table 3, we can say that the forest height estimation and the underlying ground topographic phase are reliable. Consequently, the proposed approach provides relative accuracy with small error and more accurate for vertical structural variations, especially at the sloping terrain.
To estimate the main forest parameters, the presented forest model in the alternate transmit model is used. The parameter inversion process consists of optimizing the error function and estimating the physical parameters h v , r h , f v , ϕ v ,θ, σ, f d , α, φ, ϕ d , η, f s , ϕ s , and χ. Figure 6 presents the parameter inversion performance in the 200th row of mountain forest A. The height sensitivity is given by the vertical wavenumber, which is about 0.16. This corresponds to 2π height ambiguity of about 40 m. In the experiments, the graphs display the value and the standard deviation of estimated parameters. This figure indicates that the total forest height is around 20 m [ Fig. 6(a)], the underlying topographic phase is around 0 rad [ Fig. 6(d)], the mean orientation angle of volume scattering contribution is around 65 deg [ Fig. 6(c)], and the degree of orientation randomness of volume scattering contribution is very low (σ ≈ 0) [ Fig. 6(b)].
The rotation angle of single-and double-bounce scattering components estimated by GMBD approach over mountain forest A is shown in Fig. 7. In this figure, it is shown that the large rotation angle for single-bounce scattering contribution in these areas (χ≈40 deg) and the rotation angle for double-bounce angle for double-bounce contribution is around 23 deg. This is caused by the dominance of surface scattering contribution in the sloping terrain and the increasing of the cross-polarization and off-diagonal terms.  To evaluate the proposed decomposition further, the derived scattering power for the proposed approach and the adaptive model-based decomposition method without orientation compensation corresponding to mountain forest A are shown in Fig. 8. Evidently, the results of the adaptive model-based decomposition show that the double-bounce scattering contribution (P d1 ) is too low, and the volume scattering (P v1 ) is relatively high. In the proposed decomposition, the double-bounce scattering power (P d2 ) is significantly enhanced. The power of the surface scattering component (P s2 ) is remarkably increased.
The scattering power contributions of the three scattering components are calculated and are also compared in Table 4. For forest area B, which corresponds to a relatively flat area, the two methods show relatively equivalent results. In this area, the volume scattering is dominant and very high: 94.6% for the adaptive decomposition and 96.6% for the proposed method. For mountain forest area A, for the adaptive decomposition, although the volume scattering component is still dominant, the percentage of dominant P v is reduced to 75.99% while the percentage of dominant single-and double-bounce contributions P s and P d are increased to 11.04% and 12.97%, respectively. For the proposed decomposition, even though the dominant volume scattering is still maximal, the percentage of dominant P v is significantly decreased to 62.98%. Otherwise, the percentage of surface scattering component P s is significantly increased to 20.45%. The percentage of double bounce is relatively increased to 16.57%. The reason lies in that the scattering mechanisms in patch A is strongly affected by the topographic variations. In addition, the topographic variations introduce changes in the penetration depth of microwave into the forest. Thereby, the occurrences of the ground-trunk double-bounce scattering are not as many as those of single-bounce scattering directly from the trunks or branches. Therefore, the reduced volume scattering power mainly changes into the single-bounce scattering. Figure 9 corresponds to the amplitude of the three scattering mechanisms contributing to the HH þ VV [ Fig. 9(a)] channel and HH − VV [ Fig. 9(b)] channel in the mountain forest area A. As shown, the amplitude of volume scattering contribution is dominant for all two channels, whereas the amplitude of surface scattering component is also remarkable for HH þ VV channel (as predicted by theories in Refs. 1 and 6). For HH − VV channel, the amplitude of doublebounce scattering component is too high, but it is not dominant due to slope of terrain.

Conclusions
The underestimation of the forest height and the scattering mechanism ambiguities over sloping terrain in the model-based decomposition and the height estimation method was discussed. The GMBD scheme, which uses all polarimetric interferometric coherence matrix elements, has been  proposed. All the model parameters and forest height can be adaptively obtained during the optimization procedure of the model inversion. In comparison to other methods, the proposed approach enables us to improve the estimation of forest height and ground topography over mountain forest area, as well as to retrieve additional parameters related to the degree of randomness, the main orientation of particles, the canopy layer depth, and power contribution of each scattering component. The GMBD approach is quite flexible and effective for analysis of more complex forest model over all types of terrains with PolInSAR images. Experimental results indicate that the forest parameters can be retrieved directly and more accurately using the proposed approach. In the future, further theoretical and experimental investigations will be done to improve the performance of the proposed approach.