## 1.

## 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.

## 2.

## General Scattering Model in Mountain Forest Areas for PolInSAR

## 2.1.

### 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, [${\mathrm{S}}_{1}$] and [${\mathrm{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}

## (1)

$${\overrightarrow{k}}_{Pi}=\frac{1}{\sqrt{2}}{[{S}_{hh}^{i}+{S}_{vv}^{i}{S}_{hh}^{i}-{S}_{vv}^{i}2{S}_{hv}^{i}]}^{\mathrm{T}},$$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)

## (2)

$$[T]=\langle \overrightarrow{k}{\overrightarrow{k}}^{*\mathrm{T}}\rangle =\left[\begin{array}{cc}{T}_{1}& \mathrm{\Omega}\\ {\mathrm{\Omega}}^{*T}& {T}_{2}\end{array}\right]\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{with}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\overrightarrow{k}=\left[\begin{array}{c}{\overrightarrow{k}}_{1}\\ {\overrightarrow{k}}_{2}\end{array}\right],$$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 vertical–horizontal 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

## (3)

$$\mathrm{tan}\text{\hspace{0.17em}}\psi =\frac{\mathrm{tan}\text{\hspace{0.17em}}\omega}{-\mathrm{tan}\text{\hspace{0.17em}}\gamma \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\theta}_{0}+\mathrm{sin}\text{\hspace{0.17em}}{\theta}_{0}},$$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}

## (4)

$$[\mathrm{\Omega}(\psi )]=[R(\psi )][\mathrm{\Omega}]{[R(\psi )]}^{*\mathrm{T}}=\left[\begin{array}{ccc}{\mathrm{\Omega}}_{11}(\psi )& {\mathrm{\Omega}}_{12}(\psi )& {\mathrm{\Omega}}_{13}(\psi )\\ {\mathrm{\Omega}}_{21}(\psi )& {\mathrm{\Omega}}_{22}(\psi )& {\mathrm{\Omega}}_{23}(\psi )\\ {\mathrm{\Omega}}_{31}(\psi )& {\mathrm{\Omega}}_{32}(\psi )& {\mathrm{\Omega}}_{33}(\psi )\end{array}\right],$$## (5)

$$[R(\psi )]=\left[\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{cos}\text{\hspace{0.17em}}2\psi & \mathrm{sin}\text{\hspace{0.17em}}2\psi \\ 0& -\mathrm{sin}\text{\hspace{0.17em}}2\psi & \mathrm{cos}\text{\hspace{0.17em}}2\psi \end{array}\right]\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{and}\phantom{\rule[-0.0ex]{1em}{0.0ex}}\psi =\frac{1}{4}\text{\hspace{0.17em}}\mathrm{arctan}\left\{\frac{-4\mathrm{Re}(\langle ({S}_{hh}-{S}_{vv}){S}_{hv}^{*}\rangle )}{-\langle {|{S}_{hh}-{S}_{vv}|}^{2}\rangle +4\langle {|{S}_{hv}|}^{2}\rangle}\right\},$$## 2.2.

### 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 ${\phi}_{hv}={\phi}_{h}-{\phi}_{v}$ and the interferometric phase related to the position in the vertical coordinate ${\varphi}_{S}$. Hence, the Pauli target vector of single-bounce contribution for both ends of the baseline is represented as

## (6)

$${\overrightarrow{k}}_{P1}=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\mathrm{cos}{\text{\hspace{0.17em}}\alpha}_{s}{e}^{j{\phi}_{1}}\\ \mathrm{sin}{\text{\hspace{0.17em}}\alpha}_{s}{e}^{j{\phi}_{2}}\\ 0\end{array}\right]\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}\text{and}\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}{\overrightarrow{k}}_{P2}=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{s}{e}^{j({\phi}_{1}+{\varphi}_{s})}\\ \mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{s}{e}^{j({\phi}_{1}+{\varphi}_{s})}\\ 0\end{array}\right],$$^{13}

## (7)

$${\alpha}_{s}=\mathrm{arccos}\left(\frac{|R{}_{\mathrm{H}}+{R}_{\mathrm{V}}|}{\sqrt{{|{R}_{\mathrm{H}}+{R}_{\mathrm{V}}|}^{2}+{|{R}_{\mathrm{H}}-{R}_{\mathrm{V}}|}^{2}}}\right).$$The phase ${\varphi}_{s}$ denotes the interferometric phase for the direct scattering from surface. The coefficients ${R}_{\mathrm{H}}$ and ${R}_{\mathrm{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

## (8)

$$[{T}_{s}]=\frac{{e}^{j{\varphi}_{s}}}{2}\left[\begin{array}{ccc}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}{\alpha}_{s}& \mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{s}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{s}{e}^{j{\phi}_{s}}& 0\\ \mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{s}{e}^{-j{\phi}_{s}}& {\mathrm{sin}}^{2}\text{\hspace{0.17em}}{\alpha}_{s}& 0\\ 0& 0& 0\end{array}\right],$$For a rugged terrain surface, the presence of topographic variations induces the local incidence angle $\theta $ to be a function of the local terrain slope

## (10)

$$\mathrm{cos}\text{\hspace{0.17em}}\theta =\frac{\mathrm{tan}\text{\hspace{0.17em}}\gamma \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\theta}_{0}+\mathrm{cos}\text{\hspace{0.17em}}{\theta}_{0}}{\sqrt{1+{\mathrm{tan}}^{2}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\gamma +{\mathrm{tan}}^{2}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\omega}},$$## (11)

$$[{T}_{s}(\chi )]=[R(\chi )][{T}_{s}]{[R(\chi )]}^{*\mathrm{T}}\phantom{\rule{0ex}{0ex}}=\frac{{e}^{j{\varphi}_{s}}}{2}\left[\begin{array}{ccc}2\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}{\alpha}_{s}& 2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\chi {e}^{j{\phi}_{s}}& -2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\chi {e}^{j{\phi}_{s}}\\ 2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\chi {e}^{-j{\phi}_{s}}& 2\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}2\chi & -{\mathrm{sin}}^{2}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}4\chi \\ -2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\chi {e}^{-j{\phi}_{s}}& -{\mathrm{sin}}^{2}\text{\hspace{0.17em}}{\alpha}_{s}\mathrm{sin}\text{\hspace{0.17em}}4\chi & 2{\text{\hspace{0.17em}}\mathrm{sin}}^{2}\text{\hspace{0.17em}}{\alpha}_{s}\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}2\chi \end{array}\right]\mathrm{.}$$Slope also modifies the apparent alpha parameter ${\alpha}_{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.

## 2.3.

### 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

## (12)

$${\overrightarrow{k}}_{P1}=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{D}{e}^{j{\phi}_{3}}\\ \mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{D}{e}^{j{\phi}_{4}}\\ 0\end{array}\right]\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}\text{and}\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}{\overrightarrow{k}}_{P2}=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\mathrm{cos}{\text{\hspace{0.17em}}\alpha}_{D}{e}^{j({\phi}_{3}+{\varphi}_{D})}\\ \mathrm{sin}{\text{\hspace{0.17em}}\alpha}_{D}{e}^{j({\phi}_{4}+{\varphi}_{D})}\\ 0\end{array}\right],$$## (13)

$${\alpha}_{D}=\mathrm{arccos}\left(\frac{|{R}_{\mathrm{GH}}{R}_{\mathrm{TH}}-{R}_{\mathrm{GV}}{R}_{\mathrm{TV}}|}{\sqrt{{|{R}_{\mathrm{GH}}{R}_{\mathrm{TH}}+{R}_{\mathrm{GV}}{R}_{\mathrm{TV}}|}^{2}+{|{R}_{\mathrm{GH}}{R}_{\mathrm{TH}}-{R}_{\mathrm{GV}}{R}_{\mathrm{TV}}|}^{2}}}\right).$$The coefficients ${R}_{\mathrm{GH}}$ and ${R}_{\mathrm{GV}}$ are the horizontal and vertical Fresnel reflection coefficients of the ground surface, respectively. Similarly, the vertical trunk surface has reflection coefficients ${R}_{\mathrm{TH}}$ and ${R}_{\mathrm{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

## (14)

$$[{T}_{D}]=\frac{{e}^{j{\varphi}_{D}}}{2}\left[\begin{array}{ccc}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}{\alpha}_{D}& \mathrm{cos}{\text{\hspace{0.17em}}\alpha}_{D}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{D}{e}^{j{\phi}_{D}}& 0\\ \mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{D}{e}^{-j{\phi}_{D}}& {\mathrm{sin}}^{2}\text{\hspace{0.17em}}{\alpha}_{D}& 0\\ 0& 0& 0\end{array}\right],$$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 single- and 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 $\eta $.

## (16)

$$[{T}_{D}(\eta )]=[R(\eta )][{T}_{D}]{[R(\eta )]}^{*T}\phantom{\rule{0ex}{0ex}}=\frac{{e}^{j{\phi}_{D}}}{2}\left[\begin{array}{ccc}2\text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}{\alpha}_{D}& 2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\eta {e}^{j{\phi}_{D}}& -2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\eta {e}^{j{\phi}_{D}}\\ 2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\eta {e}^{-j{\phi}_{D}}& 2{\text{\hspace{0.17em}}\mathrm{sin}}^{2}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}2\eta & -{\mathrm{sin}}^{2}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}4\eta \\ -2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\eta {e}^{-j{\phi}_{D}}& -{\mathrm{sin}}^{2}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}4\eta & 2\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}{\alpha}_{D}\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}2\eta \end{array}\right].$$For flat terrain surface, the surface and dihedral component are orthogonal.^{14} The orthogonal condition can be expressed as

## (17)

$${\alpha}_{s}+{\alpha}_{D}=\frac{\pi}{2}\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}\text{and}\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}{\phi}_{D}-{\phi}_{s}=\pm \pi .$$Therefore, the orthogonal condition reduces $({\alpha}_{S},{\alpha}_{D})$ to $\alpha $ and $({\phi}_{S},{\phi}_{D})$ to $\phi $ in Eqs. (11) and (16). Therefore, two equations [Eqs. (11) and (16)] are written as

## (18)

$$[{T}_{s}(\kappa )]=\frac{{e}^{j{\phi}_{s}}}{2}\left[\begin{array}{ccc}2\text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}\alpha & -2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\kappa {e}^{j\phi}& 2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\kappa {e}^{j\phi}\\ -2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\kappa {e}^{-j\phi}& 2\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}2\kappa & -{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}4\kappa \\ 2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\kappa {e}^{-j\phi}& -{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}4\kappa & 2\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}2\kappa \end{array}\right]$$## (19)

$$[{T}_{D}(\eta )]=\frac{{e}^{j{\phi}_{D}}}{2}\left[\begin{array}{ccc}2\text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}\alpha & 2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\eta {e}^{j\phi}& -2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\eta {e}^{j\phi}\\ 2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\eta {e}^{-j\phi}& 2\text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}2\eta & {\mathrm{cos}}^{2}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}4\eta \\ -2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\eta {e}^{-j\phi}& {\mathrm{cos}}^{2}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}4\eta & 2\text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}\alpha \text{\hspace{0.17em}}{\mathrm{sin}}^{2}\text{\hspace{0.17em}}2\eta \end{array}\right].$$## 2.4.

### 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 $\sigma $ can be evaluated in terms of $n$, as in Ref. 8. The standard deviation $\sigma $ 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

## (20)

$$[{T}_{v}(\overline{\theta},\sigma )]=[{T}_{a}]+p(\sigma )[{T}_{b}(2\overline{\theta})]+q(\sigma )[{T}_{c}(4\overline{\theta})],$$## (21)

$$[{T}_{a}]=\frac{1}{4}\left[\begin{array}{ccc}2& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right][{T}_{b}]=\frac{1}{4}\left[\begin{array}{ccc}0& 2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\overline{\theta}& 2\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\overline{\theta}\\ 2\text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\overline{\theta}& 0& 0\\ 2\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\overline{\theta}& 0& 0\end{array}\right][{T}_{c}]=\frac{1}{4}\left[\begin{array}{ccc}0& 0& 0\\ 0& \mathrm{cos}\text{\hspace{0.17em}}4\overline{\theta}& \mathrm{sin}\text{\hspace{0.17em}}4\overline{\theta}\\ 0& \mathrm{sin}\text{\hspace{0.17em}}4\overline{\theta}& -\mathrm{cos}\text{\hspace{0.17em}}4\overline{\theta}\end{array}\right].$$## 3.

## 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:

## (22)

$$[\mathrm{\Omega}(\psi )]={f}_{s}{e}^{j{\varphi}_{s}}[{T}_{s}(\chi )]+{f}_{D}{e}^{j{\varphi}_{D}}[{T}_{D}(\eta )]+{f}_{v}{e}^{j{\varphi}_{v}}[{T}_{v}(\overline{\theta},\sigma )]+[{T}_{\text{remainder}}],$$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 1718.–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}_{v}^{i}(\overline{\theta},\sigma )]$ approximates to the reference volume scattering coherence matrix by varying the randomness $\sigma $ and mean orientation angle $\overline{\theta}$ 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:

## (24)

$${f}_{v}^{i}=\left|\frac{{\mathrm{\Omega}}_{33}(\varphi )-{\mathrm{\Omega}}_{22}(\varphi )}{{T}_{v(33)}^{i}-{T}_{v(22)}^{i}}\right|$$## (25)

$${\phi}_{v}^{i}=\mathrm{arg}\left\{\frac{{\mathrm{\Omega}}_{33}(\varphi )-{\mathrm{\Omega}}_{22}(\varphi )}{{T}_{v(33)}^{i}-{T}_{v(22)}^{i}}\right\},$$After subtracting the resulting volume component, Eq. (22) becomes

## (26)

$$[{\mathrm{\Omega}}_{\text{remainder}}]=[\mathrm{\Omega}(\varphi )]-{f}_{v}^{i}{e}^{j{\phi}_{v}^{i}}[{T}_{v}^{i}(\overline{\theta},\sigma )]={f}_{s}^{i}{e}^{j{\phi}_{s}^{i}}[{T}_{s}^{i}(\kappa )]+{f}_{D}^{i}{e}^{j{\phi}_{D}^{i}}[{T}_{D}^{i}(\eta )]+[{T}_{\text{remainder}}^{i}].$$As can be seen, in matrix $[{\mathrm{\Omega}}_{\mathrm{remainder}}]$, there appear eight real unknowns $\{{f}_{s},{f}_{D},{\varphi}_{s},{\varphi}_{D},\alpha ,\phi ,\chi ,\eta \}$, a matrix $[{T}_{\text{remainder}}^{i}]$, and nine complex observables, since ${\mathrm{\Omega}}_{\mathrm{remainder}(\mathrm{i},\mathrm{j})}\ne {\mathrm{\Omega}}_{\mathrm{remainder}(\mathrm{j},\mathrm{i})}$. This formulation leads to a determined nonlinear equation system. Therefore, to determine the rest of the unknown parameters $\{{f}_{s},{f}_{D},{\varphi}_{s},{\varphi}_{D},\alpha ,\phi ,\chi ,\eta \}$ and $[{T}_{\text{remainder}}^{i}]$ 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 $\kappa $ and double bounce $\eta $ are all equal to $-\frac{\pi}{4}$ (Ref. 10). With each pair value ($\chi ,\eta $), the parameters $\{{f}_{s},{f}_{D},{\varphi}_{s},{\varphi}_{D},\alpha ,\phi \}$ are determined by employing an singular value decomposition of the ${\mathrm{\Omega}}_{\text{remainder}}$. The parameter set $\{{f}_{s},{f}_{D},{\varphi}_{s},{\varphi}_{D},\alpha ,\phi ,\chi ,\eta \}$ in this step is determined from condition minimum of Frobenius norm of matrix ${[{T}_{\text{remainder}}]}_{i}$. We show that the parameter set $\{{f}_{s},{f}_{D},{f}_{v},{\varphi}_{s},{\varphi}_{D},{\varphi}_{v},\alpha ,\phi ,\chi ,\eta \}$ 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)

## (28)

$$A=\frac{{\mathrm{\Omega}}_{12}(\varphi )\mathrm{sin}\text{\hspace{0.17em}}2\kappa +{\mathrm{\Omega}}_{13}(\varphi )\mathrm{cos}\text{\hspace{0.17em}}2\kappa}{\mathrm{cos}\text{\hspace{0.17em}}2\eta \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}2\kappa -\mathrm{sin}\text{\hspace{0.17em}}2\eta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}2\kappa}\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}\text{and}\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}B=\frac{{\mathrm{\Omega}}_{22}(\varphi )\mathrm{sin}\text{\hspace{0.17em}}4\kappa +2{\mathrm{\Omega}}_{23}(\varphi ){\mathrm{cos}}^{2}\text{\hspace{0.17em}}2\kappa}{{\mathrm{cos}}^{2}\text{\hspace{0.17em}}2\eta \text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}4\kappa -\mathrm{sin}\text{\hspace{0.17em}}4\eta \text{\hspace{0.17em}}{\mathrm{cos}}^{2}\text{\hspace{0.17em}}2\kappa}.$$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)

## (29)

$${h}_{v}=\frac{{\phi}_{v}-{\phi}_{s}}{{k}_{z}}=\mathrm{\Delta}\phi \frac{\lambda}{4\pi}\frac{R\text{\hspace{0.17em}}\mathrm{sin}\text{\hspace{0.17em}}\theta}{B\text{\hspace{0.17em}}\mathrm{cos}(\theta -\delta )},$$## 4.

## 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.

## 4.1.

### 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\times 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 three-stage 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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.

## Table 1

Forest height estimation for three approaches.

Parameters | True | Three-stage inversion | Adaptive decomposition | Proposed method |
---|---|---|---|---|

hv (m) | 18 | 13.4619 | 14.2123 | 17.2751 |

Φ0 (rad) | −0.148 | −0.1452 | −0.1350 | −0.1160 |

Average error (m) | 0 | 4.5381 | 3.7877 | 0.7294 |

RMSE (m) | 0 | 2.4290 | 4.6987 | 5.2801 |

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 $\gamma $ from 0 to 30 deg and the local ground azimuth slope $\omega $ 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 $(\gamma ,\omega )=(30\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg},45\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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.

## Table 2

Forest parameter estimation for difference slope terrains.

Parameters | True | γ=0 degω=6 deg | γ=0 degω=27 deg | γ=6 degω=18 deg | γ=18 degω=30 deg | γ=27 degω=45 deg | γ=30 degω=45 deg |
---|---|---|---|---|---|---|---|

hv (m) | 18 | 17.2096 | 17.0145 | 17.2751 | 17.1461 | 16.8917 | 16.1935 |

Φ0 (rad) | −0.148 | −0.1425 | −0.1371 | −0.1350 | −0.1022 | −0.1135 | −0.1131 |

Average error (m) | 0 | 0.7904 | 0.9855 | 0.7294 | 0.8539 | 1.1083 | 1.8065 |

RMSE (m) | 0 | 5.0174 | 5.2545 | 5.2801 | 4.8822 | 4.9909 | 3.6979 |

## 4.2.

### 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\times 10\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{m}}^{2}$. They consisted of full polarized data at L-band with 21.5-deg angle of incidence and composed of $\mathrm{12,816}\times 1129\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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 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\times 384$ and forest area B with $265\times 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:

## (30)

$${r}_{h}=\frac{2.7\text{\hspace{0.17em}}\mathrm{ln}{\overline{\text{\hspace{0.17em}}h}}_{v}-0.1}{{h}_{v}}.$$## Table 3

Parameter estimation for two approaches.

Parameter | Patch A | Patch B | ||
---|---|---|---|---|

Adaptive | Proposed | Adaptive | Proposed | |

hv (m) | 17.2845 | 19.3794 | 17.6490 | 19.3604 |

rh | 0.4394 | 0.4079 | 0.4335 | 0.4085 |

Φ0 (rad) | −0.0180 | −0.0161 | 0.0313 | 0.0302 |

RMSE (hv) (m) | 9.2848 | 8.0722 | 5.3888 | 5.0172 |

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}$, ${\varphi}_{v}$, $\overline{\theta}$, $\sigma $, ${f}_{d}$, $\alpha $, $\phi $, ${\varphi}_{d}$, $\eta $, ${f}_{s}$, ${\varphi}_{s}$, and $\chi $. 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\pi $ 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 ($\sigma \approx 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 ($\chi \approx 40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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.

## Table 4

Power of three scattering components (%) over two study areas.

Method | Mountain forest A | Forest B | ||||
---|---|---|---|---|---|---|

Ps | Pd | Pv | Ps | Pd | Pv | |

Adaptive | 11.04 | 12.97 | 75.99 | 2.5 | 2.9 | 94.6 |

Proposed | 54.19 | 16.57 | 29.24 | 1.2 | 2.2 | 96.6 |

Figure 9 corresponds to the amplitude of the three scattering mechanisms contributing to the $\mathrm{HH}+\mathrm{VV}$ [Fig. 9(a)] channel and $\mathrm{HH}-\mathrm{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 $\mathrm{HH}+\mathrm{VV}$ channel (as predicted by theories in Refs. 1 and 6). For $\mathrm{HH}-\mathrm{VV}$ channel, the amplitude of double-bounce scattering component is too high, but it is not dominant due to slope of terrain.

## 5.

## 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.

## Acknowledgments

This work was supported by National Natural Science Foundation of China under Grant No. 61271348 and National S&T Plans in Agriculture Area in the Twelfth Five-Year Period (No. 2011BAD08B02-03).

## References

S. R. CloudeK. P. Papathanassiou, “Three-stage inversion process for polarimetric SAR interferometric,” IEE Proc. Radar Sonar Navig. 150(3), 125–134 (2003), http://dx.doi.org/10.1049/ip-rsn:20030449.1350-2395Google Scholar

H. Yamadaet al., “Polarimetric SAR interferometry for forest analysis based on the ESPRIT algorithm,” IEICE Trans. Electron. 84(12), 1917–2014 (2001).IELEEJ0916-8524Google Scholar

F. GarestierT. Le Toan, “Evaluation of forest modeling for forest height estimation using P-band PolInSAR data,” in Proc. BIOGEOSAR, IEEE, Bari, Italy (2007).Google Scholar

F. GarestierT. Le Toan, “Forest modeling for forest height inversion using single baseline InSAR/PolInSAR data,” IEEE Trans. Geosci. Remote Sens. 48(3), 1528–1539 (2010), http://dx.doi.org/10.1109/TGRS.2009.2032538.IGRSD20196-2892Google Scholar

F. GarestierT. Le Toan, “Estimation of the backscatter vertical profile of a pine forest using single baseline P-band (Pol-) InSAR data,” IEEE Trans. Geosci. Remote Sens. 48(9), 3340–3348 (2010), http://dx.doi.org/10.1109/TGRS.2010.2046669.IGRSD20196-2892Google Scholar

J. D. Ballester-BermandJ. M. Lopez-Sanchez, “Applying the Freeman-Durden decomposition concept to polarimatric SAR interferometry,” IEEE Trans. Geosci. Remote Sens. 48(1), 466–479 (2010), http://dx.doi.org/10.1109/TGRS.2009.2024304.IGRSD20196-2892Google Scholar

M. NeumannL. F. FamilA. Reigber, “Estimation of forest structure, ground and canopy layer characteristics from multibaseline polarimetric interferometric SAR dat,” IEEE Trans. Geosci. Remote Sens. 48(3), 1086–1103 (2010), http://dx.doi.org/10.1109/TGRS.2009.2031101.IGRSD20196-2892Google Scholar

M. ArriJ. VanZylY. Kim, “Adaptive model-based decomposition of polarimetric SAR covariance matrix,” IEEE Trans. Geosci. Remote Sens. 49(3), 1104–1113 (2011), http://dx.doi.org/10.1109/TGRS.2010.2076285.IGRSD20196-2892Google Scholar

J. J. VanZyl, “The effect of topography on the radar scattering from vegetated areas,” IEEE Trans. Geosci. Remote Sens. 31(1), 153–160 (1993), http://dx.doi.org/10.1109/36.210456.IGRSD20196-2892Google Scholar

Y. C. LinK. Sarabandi, “Electromagnetic scattering model for a tree trunk above a tilted ground plane,” IEEE Trans. Geosci. Remote Sens. 33(4), 687–697 (1989), http://dx.doi.org/10.1109/36.406692.Google Scholar

J. S. LeeT. L. Ainsworth, “The effect of orientation angle compensation on coherence matrix and polarimetric target decomposition,” IEEE Trans. Geosci. Remote Sens. 49(1), 53–64 (2011), http://dx.doi.org/10.1109/TGRS.2010.2048333.IGRSD20196-2892Google Scholar

S. R. CloudeK.P. Papathanassiou, “Polarimetric SAR interferometry,” IEEE Trans. Geosci. Remote Sens. 36(5), 1551–1556 (1998), http://dx.doi.org/10.1109/36.718859.IGRSD20196-2892Google Scholar

I. HajnsekE. PottierS.R. Cloude, “Inversion of surface parameters from polarimetric SAR,” IEEE Trans. Geosci. Remote Sens. 41(4), 727–744 (2003), http://dx.doi.org/10.1109/TGRS.2003.810702.IGRSD20196-2892Google Scholar

G. Singhet al., “Hybrid Freeman/Eigenvalue decomposition method with extended volume scattering model,” IEEE Geosci. Remote Sens. Lett. 10(1), 81–85 (2013), http://dx.doi.org/10.1109/LGRS.2012.2193373.IGRSBY1545-598XGoogle Scholar

A. FreemanS.L. Durden, “A three component scattering model to describle polarimetric SAR data,” IEEE Trans. Geosci. Remote Sens. 36(3), 963–973 (1998), http://dx.doi.org/10.1109/36.673687.IGRSD20196-2892Google Scholar

M. ArriJ. VanZylY. Kim, “A general characterization for polarimetric scattering from vegetation canopies,” IEEE Trans. Geosci. Remote Sens. 48(9), 3349–3357 (2010), http://dx.doi.org/10.1109/TGRS.2010.2046331.IGRSD20196-2892Google Scholar

Y. Yamaguchiet al., “Four component scattering model for polarimetric SAR image decomposition,” IEEE Trans. Geosci. Remote Sens. 43(8), 1699–1706 (2005), http://dx.doi.org/10.1109/TGRS.2005.852084.IGRSD20196-2892Google Scholar

B. Zouet al., “Building parameters extraction from spaceborne PolInSAR image using a built-up area scattering model,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 6(1), 162–170 (2013), http://dx.doi.org/10.1109/JSTARS.2012.2210033.1939-1404Google Scholar

W. T. AnY. CuiJ. Yang, “Three-component model-based decomposition for polarimetric SAR data,” IEEE Trans. Geosci. Remote Sens. 48(6), 2732–2739 (2010), http://dx.doi.org/10.1109/TGRS.2010.2041242.IGRSD20196-2892Google Scholar

P.E. GillW. Murray, “Algorithms for the solution of the nonlinear least-square problems,” SIAM J. Numer. Anal. 15(5), 977–992 (1978), http://dx.doi.org/10.1137/0715063.SJNAEQ0036-1429Google Scholar

M. L. Williams, “PolSARproSim: a coherent, polarimetric SAR simulation of forest for PolSARPro,” http//earth.eo.esa.int/polsarpro/SimulatedDataSources.html (2006).Google Scholar

W. Liet al., “Impact of topography for forest height estimation,” in Proc. 1st Int. Conf. Agro-Geoinformatics, pp. 1–6, IEEE, Shanghai, China (2012).Google Scholar

K. P. PapathanassiouS. R. Cloude, “The effect of temporal decorrelation on the inversion of forest parameters from PolInSAR data,” in Proc. IEEE Int. Conf. Geoscience and Remote Sensing Symposium, pp. 1429–1431, IEEE, Toulouse, France (2003).Google Scholar

N. P. MinhB. ZouC. Yan, “Forest height estimation from PolInSAR image using adaptive decomposition method,” in Proc. 11th Int. Conf. Image Processing, pp. 1830–1834, IEEE, Beijing, China (2012).Google Scholar

## Biography

**Nghia Pham Minh** received his BS and MS degrees from the Le Qui Don Technical University, Hanoi, Vietnam, in 2005 and 2008, respectively. Now, he is pursuing his PhD degree in the School of Electronics and Information Technology at Harbin Institute of Technology (HIT), Harbin, China. He currently focuses on synthetic aperture radar (SAR) image processing, polarimetric SAR, and polarimetric SAR interferometry.

**Bin Zou** received his BS degree from the Harbin Institute of Technology (HIT), Harbin, China, in 1990, the MS degree in space studies from the International Space University, Strasbourg, France, in 1998, and the PhD degree in 2001. He is currently a professor and vice head of the Department of Information Engineering, School of Electronics and Information Technology, HIT. He currently focuses on synthetic aperture radar (SAR) image processing, polarimetric SAR, polarimetric SAR interferometry, hyperspectral imaging, and data processing.

**Hongjun Cai** received his BS degree from Zhengzhou University, Zhengzhou, China, in 2005, the MS degree from Harbin Institute of Technology (HIT), Harbin, China, in 2007, and the PhD degree from HIT, in 2012. From September 2009 to September 2010, he was a visiting student in the Satellite Geophysics Laboratory in the University of Manitoba, Canada. His research interests include polarimetric synthetic aperture radar (SAR) and polarimetric interferometric SAR data processing.

**Chengyi Wang** received his BS degree from Northeast Forestry University, Harbin, China, in 1984, and the MS degree from Department of Chinese Academy of Forestry, Beijing, China, in 1999. Recently, he has published 40 papers and participated in 20 projects. He now undertakes “Twelfth Five Year” National Science and Technology Project in Rural Areas-Research and Demonstration of Ecological Restoration Technology in the Mining Area (No. 2011BAD08B02-03).