Forest Vertical Parameter Estimation Using PolInSAR Imagery Based on Radiometric Correction

This work introduces an innovative radiometric terrain correction algorithm using PolInSAR imagery for improving forest vertical structure parameter estimation. The variance of radar backscattering caused by terrain undulation has been considered in this research by exploiting an iteration optimization procedure to improve the backscattering estimation for a Synthetic Aperture Radar (SAR) image. To eliminate the variance of backscatter coefficients caused by the local incident angle, a radiometric normalization algorithm has been investigated to compensate the influence of terrain on backscattering values, which hinders forest vertical parameter estimation. In vertical parameter estimation, species diversity and the spatial distribution of different vegetation have been modeled. Then, a combination of Fisher’s Alpha-Diversity model parameter estimation and the three-stage inversion method was designed for the vertical structure parameter. To demonstrate the efficiency of the proposed method in forest height estimation, the classical phase difference and three-stage inversion approach have been performed for the purpose of comparison. The proposed algorithm is tested on ALOS PALSAR (Advanced Land Observing Satellite Phased Array type L-band Synthetic Aperture Radar) and RADARSAT-2 (Radio Direction and Range Satellite 2) data sets for the Great Xing’an Mountain area and BioSAR (Biomass Synthetic Aperture Radar) 2007 data sets for the Remningstorp area. Height estimation results have also been validated using in-situ measurements. Experiments indicate the proposed method has the ability to compensate the influence of terrain undulation and improving the accuracy of forest vertical structure parameter estimation.


Introduction
Forest parameters, in particular vertical structure parameters, are the prerequisite of dynamic analysis of the carbon cycle and water cycle [1,2], which play an important role in research regarding environmental protection and the global climate system [3].Forest height is one of the most essential parameters of forest information.However, terrain variance poses a challenge for parameter inversion based on PolInSAR imagery because undulation may make the vegetation backscatters more complex.
Polarimetric Interferometry SAR (PolInSAR) remote sensing technology has potential in forest parameter estimation since it has sufficient penetration ability to extract the ground and canopy information of the forest.Large-scale and long time series observation of the scene can be provided to model the dynamical information of the forest area.Forest height estimation is a hot topic in PolInSAR remote sensing.
The phase difference method and model based method algorithm are the main inversion methods.The phase difference method utilizes the interferometry phase to represent the canopy and the corresponding ground for the inversion of forest height.Based on interferograms of different scattering mechanisms obtained from the coherence optimization method, Cloude [4] first uses PolInSAR technology to estimate forest height.Gabriel [5] proposes the coherence decomposition approach, where polarimetric decomposition is first used to classify different scattering mechanisms and then a phase difference procedure involving the ground phase and canopy phase is conducted to obtain forest height.Yamada [6] considers the spectral analysis theory of super-resolution technology, as well as the ESPRIT approach to estimate the orientation of signal waves, which represents the first attempt to incorporate super-resolution technology into the research on vegetation height inversion.To extract the interference phase in situations characterized by limited observation data, Wang [7] develops the ESPRIT method to calculate the interference phase and scattering power.Scattering power was utilized to classify the scattering mechanism, which uses interferograms to represent the physical characteristics of real landscape scene.Guilaso [8] suggests that ESPRIT does not exhibit accuracy in the forest terrain phase.The simulation experiment indicates that the distortion of terrain may cause unreliable estimation in the ground scattering phase, retain depolarization components and lead to Ground-over-Volume ratio of less than 3 dB.
The model based method depends on the complex coherent coefficient model of PolInSAR data, and has a strong physical and theoretical basis.In addition to the height of the vegetation, extinction coefficient and other vertical parameters can be derived.Based on the PolInSAR complex coherent coefficient model proposed by Oveisgharan [9], Papathanassiou transforms the inversion procedure into a classical non-liner parameter estimation problem, which is referred to as the six-dimension non-liner iteration method.It provides many valuable parameters simultaneously, however, calculation is time consuming, and it strongly depends on the initial value and tends to fall into the local optimal solution [10].To solve these problems, the methods such as the simulated annealing method [11], genetic algorithm [12] and BP neural network [13] are proposed.Cloude [14] proposes the three-stage inversion method under the assumption that the Ground-over-Volume ratio of different channels is located on a straight line and the estimation accuracy of the volume coherence coefficient is less than 10 dB.Thus, height estimation error can be less than 10%.Moreover, Cloude [15] proposed a multi-baseline and multi-frequency PolInSAR based method that aims at extracting the ground and canopy phase in the redundant interferometric channel to improve the accuracy of the height estimation result.However, calculating so many parameters for optimal interferometric phase extraction is time consuming.
The variance in terrain has a strong effect on vegetation height estimation.In rugged areas, terrain undulation may lead to substantial differences in ground echo, thereby changing the SAR image backscatters.As for the geometric aspect, terrain distortion changed the location of the image pixel.It may lead to shadow, layover and foreshortening in the SAR imaging area.In terms of the radiometric aspect, radar backscattering relies on the terrain since the corresponding local incident angle could change the backscattering coefficient and scattering area due to terrain variance.Thus, the land cover category will be different.As a result, even for the same land cover categories, radiometric values may differ, while different land cover categories may have a similar brightness value, which makes vegetation parameter inversion challenging work.
To eliminate distortion in the SAR imaging area, the terrain correction method is introduced via the elevation information provided by the Digital Elevation Model (DEM).The typical method includes geometric correction and radiometric correction.
Geometric correction aims at locating each pixel in the right position, in particular the real location of ridges and valleys.In terms of geometric correction methods, the Polynomial method [16] is a traditional approach that simulates the deformation area while ignoring the geometric procedure in the image area.Computation is fast but only suitable in situation of relatively flat terrain.Based on imaging theory of SAR, Leberal [17] proposes the slant distance function and Zero Doppler function for airborne SAR images.This algorithm is insufficient for space borne SAR images since the Doppler effect is not zero and distance migration is relatively significant.The DEM based simulation method first constructs the simulating image using the SAR imaging parameter and DEM information for simulation image calibration to build the coordinate function between the SAR image and DEM.Schreier [18] introduces the concept of Geocoding Ellipsoid Corrected (GEC), which transforms the orbit model, range function and Doppler frequency function into a parameterized polynomial to build the Range Doppler (RD) model.The resampling-based iteration method has also been used to obtain the real location of the SAR imagery pixel.Meier [19] introduces the Geocoded Terrain Corrected (GTC) method based on the interpolation of the satellite location and speed vector without the parameterized orbit function.In the correction step, the Doppler function-based iteration procedure is used to calculate the image coordinate from the ground coordinate without ground control points (GCP).
Radiometric correction is performed by correcting the local terrain brightness, which is defined by the incident angel, local terrain gradient and azimuth.To transform radar brightness into the backscattering coefficient, local terrain and imaging geometric should be considered.Bayer and Mauser [20][21][22] demonstrate that radar brightness can be affected by terrain in two respects.First, the scattering area may be different for the same unit in different terrain conditions, and the same ground cover may have different scattering mechanisms because of the different local incident angles induced by the different terrains.For the same scattering unit in the slant range, the forward scattering area is large, and the backward area is small in the ground range due to the terrain gradient.Small [23,24] takes full consideration of the geometric relationship between SAR imaging and DEM information in the proposed homogeneity and heterogeneity theory.As for homogeneity, the radiometric normalization method is used when the corresponding pixel in DEM has no effects on the SAR image.As for heterogeneity, the geometric relationship between the DEM and SAR image should be considered.Each DEM pixel area consists of several SAR pixels; then, the DEM pixel can be divided into several units, and the effective scattering area is represented as a combination of several units.Lars [25] proposes a project angle method for the terrain correction result in map space, which can significantly reduce the computation cost.
The different incident angles could result in huge variance in land cover backscattering.Geometric correction methods, such as DEM simulation and the RD resampling method based on interpolation, fail to achieve precision in forest parameter estimation.The radiometric correction method depends on the slope, and the terrain is relatively flat while the current method is ineffective under fractured terrain conditions.The traditional inversion procedure may limit the accuracy of vertical structure parameter estimation.Current radiometric terrain correction methods mainly focus on eliminating the influence of the amplitude of backscatter, while the phase of backscattering is ignored.
In this paper, we proposed an improved terrain correction and height estimation algorithm by considering the tree species diversity and terrain correction.Different tree species have different vertical structures and spatial distributions, and the effect of incident angle on different species is different.Land cover characteristics in SAR images are analyzed using a scattering mechanism extraction procedure via land cover classification of forest area.In terrain radiometric correction, an iteration optimization procedure is investigated for optimal backscattering estimation to compensate the influence of terrain distortion.In the vertical structure parameter extraction procedure, species diversity and the spatial distribution of vegetation are modeled.The vertical structure parameter is retrieved based on the interferometric coherence coefficient model.Compared with direct three-stage inversion, the estimation accuracy of forest height is improved.This paper is organized as follows.Section 2 describes the basic characteristics of a SAR image in forest area.A detailed framework of the proposed algorithm is presented in Section 3 after a brief introduction to the scattering mechanism classification-based radiometric correction algorithm and the vertical parameter estimation method.In Section 4, experiments to verify the proposed method and comparison methods from terrain correction analysis and tree height estimation are conducted on the real ALOS PALSAR, RADARSAT-2 and BioSAR data sets.Finally, discussions and conclusions are presented in Section 5.

PolInSAR Characteristics Analysis in Forest Area
PolInSAR image processing is a potential technique for terrain correction and vertical structure parameter extraction.In this section, PolInSAR characteristics are analyzed in terms of polarimetric SAR interferometry and backscattering.

PolInSAR Character Analysis
The PolInSAR data is acquired by two slightly different positions in space of two fully polarimetric SAR images.In the monostatic situation, each PolSAR image can be expressed as Pauli target components [26,27] Different polarization states represented by transmitting and receiving electromagnetic waves are combined to generate complex SAR data with a horizontal h or vertical v polarization basis, where i = 1, 2 indicates each one of the interferometric data sets.In the Pauli vector condition of Equation ( 1), the PolInSAR target vector is constructed from the combination of K i , for i = 1, 2.
For deterministic or point scatters, such as dihedral scattering, Equation ( 2) is a deterministic vector.As for distributed scatters, Equation ( 2) is a random vector as a consequence of the complexity of the scattering process [26].For stationary data, Equation (2) can be represented as a six-dimensional, zero-mean, complex Gaussian probability density function (pdf).Under this assumption, useful PolInSAR information can be transformed into a second-order moment as a coherency matrix.It can be reformulated as: where H indicates complex transposition and E{•} is the expectation operator.The matrices T 11 and T 22 correspond to the individual polarimetric coherency matrices, and Ω 12 is the polarimetric interferometric coherency matrix [4].
Then, the interferometric complex correlation coefficient of two complex SAR images S 1 and S 2 is obtained as Equation (4): with 0 ≤ ρ ≤ 1.In the bare region, the phase is related to the topographic height.The amplitude |ρ|, also known as the interferometric coherence, determines the accuracy of φ x .

Backscatter Conventions
In SAR imaging areas, the ratio between scattering power P s and incident power P i at ground level is defined as radar backscatter β.When the target is a distributed target, each reference area gets the backscatter ratio from the backscatter coefficient.Generally, there are three reference areas (see Figure 1).The solid rectangle A β in the slant range as the reference area and the beta naught β 0 can be expressed as [28] When the reference area is locally tangent to an ellipsoidal model of the ground surface A σ in the ground area, the result is sigma naught σ 0 E .
When the reference area is defined to be in the plane perpendicular to the line of sight from the sensor to an ellipsoidal model of the ground surface A γ , gamma naught γ 0 E is the result.
The σ 0 E or γ 0 E backscatter values can be terrain-geocoded using a digital height model (DEM), i.e., resampled into a map geometry, producing a geocoded-terrain-corrected (GTC) product [29].Although the position or geometry of the backscatter estimate has been corrected in GTC products, the radiometry of the resulting image remains ellipsoid-model based.
ISPRS Int.J. Geo-Inf.2016, 5, 186 5 of 21 When the reference area is defined to be in the plane perpendicular to the line of sight from the sensor to an ellipsoidal model of the ground surface Aγ, gamma naught   0 is the result.
The   0 .or   0 backscatter values can be terrain-geocoded using a digital height model (DEM), i.e., resampled into a map geometry, producing a geocoded-terrain-corrected (GTC) product [29].Although the position or geometry of the backscatter estimate has been corrected in GTC products, the radiometry of the resulting image remains ellipsoid-model based.

Iteration Based Backscattering Coefficient Optimization
To eliminate the effects of terrain induced geometric and radiometric distortion and maximize the correlation with the biophysical parameter as well as minimize the root mean square error of vegetation parameter estimation, a scattering mechanism classification-based radiometric terrain correction is introduced to improve the accuracy of backscattering estimation.Let polarimetric covariance matrices   ∈  × be given for i ∈ [1, … , n] , where P is the dimensionality of polarimetry.In the case of conventional polarimetric SAR measurements, the backscattering at an arbitrary combination of transmit and receive polarization is given by: where ω ∈   is a singular weighting vector  † ω and † represents the Hermitian transpose.For quad-pol polarization basis P = 4, ω can be represented by the electromagnetic wave polarization orientation and ellipticity angles associated with the sphere, Ψ and Χ for the receiving and transmitting polarizations r, t.
  , , , /,1/2 are the first and second elements of the unitary wave polarization vectors for the receive and transmit polarizations, respectively Terrain of vegetation area in SAR image: (a) SAR imaging geometry; and (b) SAR backscattering.

Iteration Based Backscattering Coefficient Optimization
To eliminate the effects of terrain induced geometric and radiometric distortion and maximize the correlation with the biophysical parameter as well as minimize the root mean square error of vegetation parameter estimation, a scattering mechanism classification-based radiometric terrain correction is introduced to improve the accuracy of backscattering estimation.Let polarimetric covariance matrices C i ∈ C P×P be given for i ∈ [1, . . . ,n], where P is the dimensionality of polarimetry.In the case of conventional polarimetric SAR measurements, the backscattering at an arbitrary combination of transmit and receive polarization is given by: where ω ∈ C P is a singular weighting vector ω † ω and † represents the Hermitian transpose.
For quad-pol polarization basis P = 4, ω can be represented by the electromagnetic wave polarization orientation and ellipticity angles associated with the sphere, Ψ and X for the receiving and transmitting polarizations r, t.
P t/r,1/2 are the first and second elements of the unitary wave polarization vectors for the receive and transmit polarizations, respectively In this case, ω can be modeled as a selection of the scattering mechanism type, including isotropy versus anisotropy, even versus odd bounce, and horizontally versus vertically oriented.Thus, the selection of polarization affects the sensitivity to a specific scattering mechanism in the illuminated resolution cell, which has a strong correlation with the terrain.The polarization basis of C i can vary simultaneously by applying the singular transformation matrix U ∈ C P × P , which can be rewritten as: To eliminate the distortion induced by terrain, the key problem is to find the polarization state that is the most sensitive to the effect of vegetation, which can be formulated by: (1) maximizing the absolute value of correlation (R) between the backscattering and the prior polarization information with the biophysical parameter from the tie points of the original measurements; (2) maximizing the coefficient of sensitivity factor between the specific scattering mechanism and backscattering (R 2 ); and (3) minimizing the root mean squared error of the inverted estimates (RMSE).If the relationship between the backscattering coefficient and polarimetric from the prior information is linear, the optimal procedures to maximize R and R 2 and minimize RMSE for estimating the weight parameter ω are unbiased, and all three optimal problems are equivalent.The correlation between polarization information y i , which indicates the relationship between land cover and corresponding backscattering in the SAR image, and the backscatter coefficient σ i(ω) is described as Thus, the polarization-dependent optimization of the absolute value of correlation can be formulated as where the arg max () operation stands for the argument of the maximum and A are Hermitian positive semi-definite, while B i is Hermitian.To solve the optimization problem of Equation ( 13), an iterative procedure is used to obtain the approximate solution.The principle of the optimization procedure is to transform C i to the optimal polarization and optimize an approximate correlation function in the new polarization basis, using inequality Equation ( 14): Then, the optimization function in Equation ( 13) can be represented as below: Equation ( 15) is used to maximize the lower boundary of the function, where F is a unique Hermitian positive semi-definite square root of ∑ B † i B i .To simplify the calculation, we neglect the absolute value operation in Equation (15).Then Lagrangian multipliers can be used to solve the optimal problem as the following: We obtain the optimization result by maximizing the numerator and keeping the denominator as a constant of Equation ( 15) simultaneously.If the sensitivity ω † Aω is negative, it may lead to the optimal correlation being less than zero.In this case, the matrix A should be multiplied by −1 to find the optimal polarization.With respect to the derivation of L in Equation ( 16), the solution can be transformed into an Eigen-problem.
The problem formulation and solution in Equations ( 16) and ( 17) provide an initial approximation, which is close to the optimum.To improve the result further, an iterative approach is investigated.It is based on alternatively optimizing L in Equation ( 16) and transforming the polarization basis of C i toward a state with the optimal polarization, as shown in Equation ( 8).As a result, we can optimize the estimation of the backscattering coefficient by solving problems that are sensitive to terrain induced distortion.

Radiometric Based Terrain Correction
Terrain undulation induced by the local incident angle in PolInSAR imagery leads to a decrease in the interferometric, which may result in an unstable estimation of height in forest area.To solve this problem, radiometric based terrain correction is proposed.The procedure of the proposed method is demonstrated using the following steps: 1. Pre-processing: A radiometric calibration step is used to convert the original SAR image into radar brightness β, and the multi-look and filtering process is conducted to obtain the slant SAR data.2. DEM area calculation: Due to the fact that the edge points of DEM have different elevations, the projection area on the ground range corresponding to the DEM grid can be represented by the summation of two triangles in Figure 2a for accurate projection area calculation.In the experiment, we choose ASTER DEM with a resolution of 30 m, downloaded from http://gdem.ersdac.jspacesystems.or.jp/. ) and projection calculation on the ground range: First, this procedure calculates the original SAR image coordinate (r init , a init ) and local incident angle corresponding to the DEM cell.Then, the projections on the γ plane dγ DEM (N, E) are achieved based on the local incident angle.Finally, the SAR coordinate and local incident angle are combined to build the LUT. 4. SAR image simulation: The original simulation result of the SAR image can be obtained by using the LUT and the projections on the γ plane.Since the coordinates in the LUT were not integers, the area should be assigned to the neighbor pixel cell.Then, it can be realized through a bilinear interpolation procedure shown in Figure 2b, where (int) means integer operation.

LUT (Look Up Table
5. Registration of simulation and real SAR image: We first calculate the offset between the original SAR image and simulation SAR image to build a registration polynomial.Then, a cross searching procedure is performed to obtain the match point.Next, the match point is used to build the quadratic polynomial function to revise the LUT.Finally, we can get the modified simulation SAR image from the γ projection and revised LUT.

Radiometric correction:
The radiometric correction result on the slant range can be obtained by the area integral on the γ plane.Then, the Orthophoto of the final terrain correction result is obtained by using the revised LUT through a bilinear interpolation procedure.
5. Registration of simulation and real SAR image: We first calculate the offset between the original SAR image and simulation SAR image to build a registration polynomial.Then, a cross searching procedure is performed to obtain the match point.Next, the match point is used to build the quadratic polynomial function to revise the LUT.Finally, we can get the modified simulation SAR image from the γ projection and revised LUT. 6. Radiometric correction: The radiometric correction result on the slant range can be obtained by the area integral on the γ plane.Then, the Orthophoto of the final terrain correction result is obtained by using the revised LUT through a bilinear interpolation procedure.

RVoG Model Based Forest Vertical Structure Modeling
For PolInSAR based vegetation height estimation, the phase difference method failed to extract the real phase center of the canopy and ground; it tends to underestimate the real height.The six-dimension-based non-liner iteration approach has a strong dependence on the original parameter and the selection of the polarization channel.Moreover, it is often time consuming and lacks a physical mechanism.Thus, we introduced the improved RVoG model based three-stage inversion procedure to estimate the vegetation height.
The most important aspect of PolInSAR data analysis is the dependency of the interferometric information on polarization.To derive the underlying topography, the polarimetric effects of the ground, including the direct and the double-bounce effects, rely on the principle that the doublebounce effect of the polarimetric phase has opposite signs in different terms of the polarimetric interferometric covariance matrix, while the interferometric phase appears with the same sign [26].
Forest scattering is a two-layer model, shown in Figure 3.One is volume scattering from the canopy and the other is surface scattering from the ground.To acquire a better estimation of the scattering center of the ground and volume scattering, the scattering mechanism classification method is exploited.For volume scattering, it can be separated if the projection vector of the two fully polarimetric images is properly selected.

RVoG Model Based Forest Vertical Structure Modeling
For PolInSAR based vegetation height estimation, the phase difference method failed to extract the real phase center of the canopy and ground; it tends to underestimate the real height.The six-dimension-based non-liner iteration approach has a strong dependence on the original parameter and the selection of the polarization channel.Moreover, it is often time consuming and lacks a physical mechanism.Thus, we introduced the improved RVoG model based three-stage inversion procedure to estimate the vegetation height.
The most important aspect of PolInSAR data analysis is the dependency of the interferometric information on polarization.To derive the underlying topography, the polarimetric effects of the ground, including the direct and the double-bounce effects, rely on the principle that the double-bounce effect of the polarimetric phase has opposite signs in different terms of the polarimetric interferometric covariance matrix, while the interferometric phase appears with the same sign [26].
Forest scattering is a two-layer model, shown in Figure 3.One is volume scattering from the canopy and the other is surface scattering from the ground.To acquire a better estimation of the scattering center of the ground and volume scattering, the scattering mechanism classification method is exploited.For volume scattering, it can be separated if the projection vector of the two fully polarimetric images is properly selected.
where ω V 1 and ω V 2 mean the projection vector and T V 2 represents the volume scattering of the image.To separate the surface scattering Tg , the extended Freemen-Durden [30] decomposition is used.
where T 11 (i, j) represent the element of column i and row j of the polarimetric coherency matrix T 11 .
where T11(i, j) represent the element of column i and row j of the polarimetric coherency matrix T11.After the scattering mechanism classification procedure, the complex interferometric coherence for a random volume over ground (RVoG) is derived as shown in Equations ( 23) and (24), which is a three-component singular complex vector defining the choice of polarization [4], σ is the extinction coefficient of the medium, kz is the vertical wave number of the interferometer (following spectral range filtering) and θ is the incident angle.The angles φ1 and φ2 are the phase centers of the bottom of layers 1 and 2, respectively.
Equation (24) indicates that the complex coherence coefficient in the complex plane should be located in a straight line; the two endpoints are γ ∅ 0 and  ∅ 0 respectively.Since m(ω) equals zero, the pure vegetation volume coherence coefficient is obtained, and it can be used to estimate tree height and the extinction coefficient in three-stage inversion.After the scattering mechanism classification procedure, the complex interferometric coherence for a random volume over ground (RVoG) is derived as shown in Equations ( 23) and ( 24), which is a three-component singular complex vector defining the choice of polarization [4], σ is the extinction coefficient of the medium, k z is the vertical wave number of the interferometer (following spectral range filtering) and θ is the incident angle.The angles ϕ 1 and ϕ 2 are the phase centers of the bottom of layers 1 and 2, respectively.

Framework of Terrain Correction and Species Diversity Based Vegetation Vertical Structure Inversion Algorithm
Equation (24) indicates that the complex coherence coefficient in the complex plane should be located in a straight line; the two endpoints are γe i∅ 0 and e i∅ 0 respectively.Since m(ω) equals zero, the pure vegetation volume coherence coefficient is obtained, and it can be used to estimate tree height and the extinction coefficient in three-stage inversion.

Framework of Terrain Correction and Species Diversity Based Vegetation Vertical Structure Inversion Algorithm
The framework of the proposed method is shown in Figure 4. First, the scattering mechanism classification procedure is exploited based on the PolInSAR data for land cover scattering classification.Then, the iteration-based backscattering optimal method is investigated to obtain the backscattering coefficient which is the most sensitive to terrain distortion.The details are described in Section 3.1.Then, the radiometric terrain correction method is investigated to eliminate the effect of terrain undulation described in Section 3.2.Finally, after Fisher's alpha-Diversity estimation, the RVoG model based three-stage inversion procedure is introduced to obtain the vegetation vertical structure parameters.classification.Then, the iteration-based backscattering optimal method is investigated to obtain the backscattering coefficient which is the most sensitive to terrain distortion.The details are described in Section 3.1.Then, the radiometric terrain correction method is investigated to eliminate the effect of terrain undulation described in Section 3.2.Finally, after Fisher's alpha-Diversity estimation, the RVoG model based three-stage inversion procedure is introduced to obtain the vegetation vertical structure parameters.Based on the distribution character of the RVoG complex coherence model, the main procedure of the three-stage inversion algorithm contains three steps: 1. Least squares straight line fitting: Based on the complex coherence coefficient of different polarization channels, the least squares method is used to fit the straight line.If there are only two polarization channels, fitting process is degenerated for calculating the straight line of the two points.2. Terrain phase calculation: the corresponding complex coherence coefficient of the terrain interferometry phase depends on the intersection of the unit circle and fitting line.It can be obtained by calculating the relative location between the channel with the largest volume scattering and the channel with the largest surface scattering.After the optimal of the backscatter coefficient, the terrain correction result is calculated with gamma naught and the backscattering area A. 3. Vegetation height estimation: Based on the terrain phase in Step 2, the multiplied complex coherence coefficient γ with e −iφ eliminates the terrain phase to obtain pure volume scattering coherence coefficient γV.Then, the vegetation height can be calculated as below: Due to the variance of distribution and the structure of different forms of vegetation, the accuracy of the retrieval result is restricted.A diversity coefficient based on Fisher's Alpha-Diversity decomposition [31] is proposed for the vertical structure parameter, which can be calculated individually among different species.The alpha factor can be used to describe the relationship between Based on the distribution character of the RVoG complex coherence model, the main procedure of the three-stage inversion algorithm contains three steps: 1.
Least squares straight line fitting: Based on the complex coherence coefficient of different polarization channels, the least squares method is used to fit the straight line.If there are only two polarization channels, fitting process is degenerated for calculating the straight line of the two points.

2.
Terrain phase calculation: the corresponding complex coherence coefficient of the terrain interferometry phase depends on the intersection of the unit circle and fitting line.It can be obtained by calculating the relative location between the channel with the largest volume scattering and the channel with the largest surface scattering.After the optimal of the backscatter coefficient, the terrain correction result is calculated with gamma naught and the backscattering area A.

3.
Vegetation height estimation: Based on the terrain phase in Step 2, the multiplied complex coherence coefficient γ with e −iϕ eliminates the terrain phase to obtain pure volume scattering coherence coefficient γ V .Then, the vegetation height can be calculated as below: Due to the variance of distribution and the structure of different forms of vegetation, the accuracy of the retrieval result is restricted.A diversity coefficient based on Fisher's Alpha-Diversity decomposition [31] is proposed for the vertical structure parameter, which can be calculated individually among different species.The alpha factor can be used to describe the relationship between the number of species and the number of individuals in the image scene, which indicates the variance in spatial distribution and attribution of different species as in Equation ( 26): where S is the number of species.N is the number of individuals, and α is the diversity coefficient.
The alpha factor can be used to represent the spatial distribution and forest vegetation diversity of the forest area using in-situ measurement results.Then, the estimation of vegetation height is transformed into an optimization problem as below:

Experimental Data and Description
In this section, radiometric terrain correction and the height estimation algorithm are tested on three data sets that include ALOS PALSAR, RADARSAT-2 and BioSAR 2007 imagery.Detailed information regarding image size and location are shown in Table 1.For ALOS PALSAR quad polarization SAR data, the spatial resolution is 30 m.It was acquired in April 2011 in the Great Xing'an Mountain area in Genhe, China, and its image size is 1248 × 18,432 pixels.The original data set and corresponding optical image from Google earth are presented in Figure 5a,b, respectively.RADARSAT-2 data were also acquired in quad polarization mode with a resolution of 4.7 m × 5.1 m from May to August 2013, in the Great Xing'an Mountain Yigen, China.Its image size is 3572 × 5914 pixels.The main land cover is grass and crop.The average altitude of the Great Xing'an Mountain area is over 1000 m with a slope angle from 20 to 50 degrees.This area has complex terrain characteristics of undulation and fragments; the relative altitude of this area ranges from 100 to 300 m.In addition, the main ground cover categories are forests and grasses with a fraction rate of nearly 75%.The dominant species are larix and white birch.BioSAR 2007 PolInSAR data with P and L bands were acquired from March to April 2007 with the DLR E-SAR sensor over the Remningstorp area in the south area of Sweden.The spatial resolution is 2 m × 2 m.The altitude of this forest area ranges from 120 to 145 m.Its image size is 13,641 × 1483 pixels, and the dominant species are spruce, pine and birch.The average height of the vegetation in different test locations ranges from 10 m to 30 m.This area has been extensively studied and in-situ measurements were conducted, including, the size of the area, vegetation height, fraction, and the main species, during acquisitions in 17 areas in the image scene to evaluate the results of different inversion methods, as shown in Table 2.For terrain correction, mean value, standard deviation, slope between local incident angle and backscattering, correction rate and calculation time are used to evaluate the performance of different terrain correction algorithm.The correction rate Cσ is defined as the decreasing percentage of the standard deviation between the original SAR imagery and terrain corrected SAR imagery, as seen in Equation (28).Sa and Sb represent the standard deviation of the original SAR image and terrain corrected image, respectively.For vegetation height estimation, mean value and standard deviation are also used to evaluate the performance of different algorithms.The proposed method is run in the MATLAB 2011a platform and polarimetric decomposition and interferometric are run in PolSARpro software.Comparison methods including terrain correction are conducted using Gamma software with an Intel Core CPU @2.4 GHz and 48.00 GB RAM.For terrain correction, mean value, standard deviation, slope between local incident angle and backscattering, correction rate and calculation time are used to evaluate the performance of different terrain correction algorithm.The correction rate C σ is defined as the decreasing percentage of the standard deviation between the original SAR imagery and terrain corrected SAR imagery, as seen in Equation (28).S a and S b represent the standard deviation of the original SAR image and terrain corrected image, respectively.For vegetation height estimation, mean value and standard deviation are also used to evaluate the performance of different algorithms.The proposed method is run in the MATLAB 2011a platform and polarimetric decomposition and interferometric are run in PolSARpro software.Comparison methods including terrain correction are conducted using Gamma software with an Intel Core CPU @2.4 GHz and 48.00 GB RAM.

Terrain Correction Result
Table 3, Figures 6 and 7 presents the results of different terrain correction algorithms as applied to the ALOS PALSAR data for Genhe, China.To compare the performance of the test algorithms, the scatter plots in Figure 7 are presented to show the slope between the incident angle and backscatter coefficient of different terrain correction methods.In total, 100,000 points are selected using the Gaussian random sampling method from the original SAR image.Then, the least squares method is investigated to obtain the slope between the local incident angle and backscatter coefficient.To eliminate the random error of sample selection, these procedures are repeated five times, and the final slope is the average value.From the original SAR data in Figure 6a, we can see that the mean value of backscattering is the lowest, with the largest standard deviation of 4.9883 dB, which indicates that the effects of terrain undulation decrease the backscattering signal and increase the uncertainty of backscattering estimation.The terrain undulation and backscattering coefficient is strongly correlated since the backscattering value is decreasing with the increase in the local incident angle in the original SAR imagery, and the slope is −0.1625.For the RD method, although the mean value of backscattering is increasing when compared with the original SAR image, the standard deviation is still over 4 dB; we can see some "high light" area in Figure 6b.Unlike the RD method, the radiometric normalization approach restrains the "light" phenomenon in the SAR image; the slope is −0.0845.Thus, the correction rate is increasing, but the performance has a slight increase in the slope area, which suggests that local terrain and backscattering are still correlated.The proposed method considers the effects of terrain induced variance of the incident angle, so the correlation between the local incident angle and backscattering diminishes.As a result, it has the largest slope of −0.0322, and the correction rate on the proposed method is over 30%.The computation time of the RD based method, radiometric normalization method and the proposed method on the ALOS PALSAR data set is 23, 28 and 35 min, respectively.To investigate the statistical significance of the sample selection, we calculated the slope using different sample points including 1000, 10,000, 100,000 and 1,000,000 and whole points of the imagery in the ALOS PALSAR data with the four methods.The mean slope and relative error are illustrated in Table 4.The mean slope is the mean value of the slope using different sample points in five instances, and the relative error is calculated as in Equation ( 29) below, where s sample means the mean value of the slope using different sample points and s img represents the slope calculated using the whole points of the image.The results in Table 4 indicate that the accuracy of the slope is increasing with the increase in sample points.The relative error using 100,000 and 1,000,000 is less than 5%, which indicates a higher statistical significance with respect to the whole data set.More sample point means a greater calculation time, so we choose 100,000 sample points to calculate the slope in our experiments.
The results of the proposed method on the RADARSAT-2 image with four different acquisition times are illustrated in Figures 8 and 9 with a calculation time of 32 min on each scene.The scatter plots in Figure 9 are obtained using the same method as in Figure 7. Table 5 gives the quantative results of terrain correction on the proposed algorithm.We can see that all correction rates on the four images are greater than 20%, with a slope of around −0.06.The result for data obtained on 23 May 2013 had the largest slope when compared to the other three images, which means that the backscattering coefficients are less correlated to the local incident angle.There are many factors that may cause the variance in the backscattering coefficients.June to August 2013 was the growth period for grasses and crops, which may explain why the slope in Figure 9a is larger than the slope of the other three results.Soil moisture is another factor leading to the variance in backscattering.The meteorological data reported that it did not rain during the SAR data acquisition time, so the soil moisture may not be the main factor that causes the variance in backscattering.This experiment indicates that vegetation parameters are also correlated to backscattering, which may also influence the result of terrain correction.Since the orbit of the RADARSAT-2 sensor is not stable, we could not obtain the stable interferometric, and this represents a lack of physical support when used for height estimation.For the ALOS PALSAR data set, we did not obtain the in-situ measurements during the acquisition time.For BioSAR 2007 data, which provide the in-situ measurements, we only carried out height estimation analysis in BioSAR data with the P and L band.Forest height estimation experiments were conducted on this data set to demonstrate the performance of different methods and the performance of the L and P band PolInSAR data.

Vegetation Height Estimation Result
Since the orbit of the RADARSAT-2 sensor is not stable, we could not obtain the stable interferometric, and this represents a lack of physical support when used for height estimation.For the ALOS PALSAR data set, we did not obtain the in-situ measurements during the acquisition time.For BioSAR 2007 data, which provide the in-situ measurements, we only carried out height estimation analysis in BioSAR data with the P and L band.Forest height estimation experiments were conducted on this data set to demonstrate the performance of different methods and the performance of the L and P band PolInSAR data.Figure 10 shows the results of vegetation height estimation for the BioSAR data set.The height is divided into three intervals with an average relative height of 6-12 m, 12.1-18 m and 18.1-24 m based on the in-situ measurement.Then, the height estimation results based on phase difference, three-stage inversion and the proposed method are illustrated in Figure 10 and Table 6.The phase difference method used the lowest time of 0.8 h, followed by the three-stage inversion method with 1.1 h.However, the proposed method requires more time to calculate the species factor with a computation time of 1.2 h.We can see that all three approaches underestimate height, and the standard deviation increases with the increase in forest height.The phase difference method obtained the lowest height estimation result of the three methods since the phase center corresponding to the ground and canopy scattering mechanism was inaccurate.It is usually characterized by underestimation of the canopy phase center and overestimation of the ground phase center.The three-stage inversion method considers the physical mechanism between the ground and canopy in the PolInSAR data; the height estimation results improved when compared with those of the phase difference method.The comparison between the three-stage inversion method and the proposed method indicates that the variance in spatial distribution and attribution variance can affect the result of height estimation.For the proposed method, the mean value of the height estimation has the lowest bias when compared with in-situ measurements, and the standard deviations of the three methods are nearly the same.Moreover, we also compared L band and P band height estimation results using the proposed method, as shown in Figure 11 and Table 7.It is clear that the bias of the height estimation is increasing with the increase in tree height in both the P and L band data, and the calculation time is 1.2 h in both images.In 6-12 m and 12.1-18 m range, the standard deviation is just around 1 m.For 18.1-24 m, the standard deviation is nearly 3 m.This may be caused by the calculation error of the coherence coefficient induced by the penetration depth and extinction coefficient.Since the penetration depth of the P band is larger than that of the L band, we can obtain a better estimation of the terrain and canopy scattering phase center, and the height estimation results for the P band are better than those of the L band with a higher mean value and lower standard deviation.

Conclusions
This paper presents a radiometric terrain correction method and forest vertical structure parameter estimation approach in PolInSAR imagery.Terrain undulation is considered, and a form of scattering mechanism classification and a scattering center optimization method have been proposed to improve the accuracy of backscattering coefficient estimation.The purpose is to eliminate the influence of complex terrain.In the forest height estimation procedure, the species diversity factor is exploited in the inversion model to describe the spatial distribution as well as the scattering variance of different species.To validate the effectiveness of the scattering center optimization approach for terrain correction, a Geometric-based method and radiometric normalization method have been applied to the on ALOS PALSAR and RADARSAT-2 data sets.The result shows that the proposed method obtains the highest mean value and lowest standard deviation of the backscattering coefficient, and the correction rate of the proposed method is also the highest.In the tree height estimation experiment, the proposed method obtain get the most accurate result for tree height at different height intervals, and the P band inversion result outperforms the L band inversion result with higher height estimation and lower standard deviation.The main contributions of the proposed algorithm include two aspects.A scattering mechanism classification and scattering center optimization approach is proposed to eliminate the influence of terrain undulation.The spatial distribution and scattering variance of different species have been considered in the height inversion procedure to improve the accuracy of height estimation.Further work will be extended to analyze the effect of surface soil moisture on backscattering and carbon cycling will be modeled in the forest area, which has been hindered by terrain undulation.A robust and operational model will be constructed for analyzing the relationship between biomass estimation and carbon cycling in the forest area.

Figure 1 .
Figure 1.Terrain of vegetation area in SAR image: (a) SAR imaging geometry; and (b) SAR backscattering.

Figure 3 .
Figure 3. Two-layer coherence RVoG model illustration for vegetated land surface.The different colors of the arrow on the left represent different scattering mechanisms and the corresponding phases are shown on the right side.The surface, double-bounce, and short understory cause similar interferometric phases at the ground level.The volume layer introduces a phase offset and an additional decorrelation due to phase variation.

Figure 3 .
Figure 3. Two-layer coherence RVoG model illustration for vegetated land surface.The different colors of the arrow on the left represent different scattering mechanisms and the corresponding phases are shown on the right side.The surface, double-bounce, and short understory cause similar interferometric phases at the ground level.The volume layer introduces a phase offset and an additional decorrelation due to phase variation.

Figure 4 .
Figure 4. Framework of terrain correction and species diversity based vegetation vertical structure inversion.

Figure 4 .
Figure 4. Framework of terrain correction and species diversity based vegetation vertical structure inversion.

Figure 5 .
Figure 5. Test data sets in Genhe, the Great Xing'an Mountain area, China: (a) ALOS PALSAR data; and (b) optical data.

Figure 5 .
Figure 5. Test data sets in Genhe, the Great Xing'an Mountain area, China: (a) ALOS PALSAR data; and (b) optical data.

Figure 7 .
Figure 7.Comparison of backscattering coefficient between measurement data and experiment result on the ALOS PALSAR data set for Genhe: (a) experiment data; (b) RD based geometric correction; (c) radiometric normalization correction; and (d) backscatter optimization radiometric correction.

Figure 11 .
Figure 11.Height estimation of the proposed method on the BioSAR 2007 with P and L band for Remningstorp: (a) P band; and (b) L band.

Table 1 .
Overview of experimental SAR imagery.

Table 2 .
Description of the in-situ measurements of forests in the Remningstorp area.

Table 2 .
Description of the in-situ measurements of forests in the Remningstorp area.

Table 3 .
Backscatter coefficient comparison of different terrain correction method on the ALOS PALSAR data set in Genhe.

Table 4 .
The statistical significance of sample selection.

Table 3 .
Backscatter coefficient comparison of different terrain correction method on the ALOS PALSAR data set in Genhe.

Table 4 .
The statistical significance of sample selection.

Table 5 .
Time series analysis of backscatter coefficient on RADARSAT-2 data sets for Yigen.Calculation time: 32 min.

Table 5 .
Time series analysis of backscatter coefficient on RADARSAT-2 data sets for Yigen.Calculation time: 32 min.

Table 5 .
Time series analysis of backscatter coefficient on RADARSAT-2 data sets for Yigen.Calculation time: 32 min.

Table 6 .
Accuracy results of different height estimation methods on BioSAR 2007 data set for Remningstorp.

Table 7 .
Accuracy results of the proposed method with P and L band BioSAR data set for Remningstorp.