S-RVoG Model Inversion Based on Time-Frequency Optimization for P-Band Polarimetric SAR Interferometry

Abstract

This paper investigates the potential of the time-frequency optimization on the basis of the sublook decomposition for forest height estimation. The optimization is deemed to be capable of extracting a relatively accurate volume contribution when P-band polarimetric interferometric synthetic aperture radar (Pol-InSAR) systems are adopted to observe forest-covered areas. The highest and the lowest phase centers acquired by the time-frequency optimization modify the conventional three-stage inversion process. This paper presents, for the first time, a performance assessment of the time-frequency optimization on P-band Pol-InSAR data over boreal forests. Simultaneously, to alleviate the model inversion errors caused by topographic fluctuations, forest height is estimated based on the sloped Random Volume over Ground (S-RVoG) model in which the incidence angle is corrected with the terrain slope. The E-SAR P-band Pol-InSAR data acquired during the BIOSAR 2008 campaign in Northern Sweden is utilized to evaluate the performance of the proposed method. From the results of the forest height estimation preprocessed with time-frequency optimization, the root mean square error (RMSE) of Random Volume over Ground (RVoG) and S-RVoG model on negative slope are 5.09 m and 4.71 m, respectively. It is concluded that the time-frequency processing and negative terrain slope compensation improve the inversion performance by $41.49\%$ and $11.96\%$, respectively.

1. Introduction

On account of the significant penetration of P-band polarimetric interferometric synthetic aperture radar (Pol-InSAR) systems, all polarimetric channels contain a non-ignorable ground contribution [1,2,3,4,5], generally leading to an underestimation of forest height in view of the Random Volume over Ground (RVoG) model [3,6]. So far, the RVoG model depicting the forest vertical backscatter profile with an exponential function has been widely taken advantage of in forest height estimation [7,8,9,10,11,12,13,14,15,16,17,18,19,20]. Moroever, when referring to temporal decorrelation [21,22,23,24,25], canopy filling [12,26] and terrain fluctuations [6,27,28,29,30,31,32], some modified models were proposed to complement the RVoG model. Under the circumstance of the Quad-pol single-baseline observation, the forest height estimation based on the RVoG model is a six-dimensional parameter optimization problem [1,3,5,10,12]. In order to retrieve the forest height with Quad-pol single-baseline data, the null ground-to-volume ratio (NGVR) hypothesis is frequently adopted to derive a unique solution from the RVoG inversion equation [3,12]. Specifically, the NGVR hypothesis indicates that the ground-to-volume ratio (GVR) of at least one polarimetric channel is sufficiently small, usually less than −10 dB [12]. In view of the RVoG model, the three-stage [12] inversion process applying single-baseline data has been widely implemented for forest height acquisition [15,31,33,34]. That is a geometric inversion method that requires fitting the estimated coherences of various polarimetric channels to obtain the coherence line on the complex plane. Thus, the ground level and canopy contributions deemed to be on the coherence line can then be achieved. Ultimately, the forest height is extracted from the estimated volume coherence by the look-up table method.

Previous work has manifested that the forest height estimation based on RVoG model presents a relatively high accuracy at L-band Pol-InSAR observations over boreal forests [6]. However, due to the strong penetration of P-band and the sparsity of boreal forest, the traditional three-stage method based on the RVoG model is no longer effective, specifically, leading to a serious underestimation. That is because the NGVR assumption is invalid under this circumstance. Up to now, two main strategies have been adopted to ameliorate the underestimation [3,35]. The first strategy is to increase the observation space with multi-baseline data [3,26,36,37,38]. Nevertheless, it brings about a high complexity of data acquisition and model inversion. The second strategy is to fix the extinction coefficient so that the single-baseline inversion is still feasible [1], yet the inversion performance is heavily dependent on the choice of the fixed extinction coefficient. To overcome the above limitations, the time-frequency optimization was developed as a method independent of the above two strategies, allowing a separation of the canopy and ground contributions with single-baseline Pol-InSAR data [2]. The sublook decomposition is operated in the azimuth direction on the Pauli and lexicographic base polarimetric channels. Then, the polarimetric channels corresponding to the highest and the lowest phase centers are assumed to be related to the canopy and ground contributions, respectively. Since it exclusively takes into account the polarimetric channels consistent with the highest and lowest phase centers, with different polarimetric coherences not necessarily located on a straight line on the complex plane, this method fundamentally derives from the Oriented Volume over Ground (OVoG) [9,39,40,41] model. The ground coherence is not constrained to equal one, as the lowest phase estimated by time-frequency optimization is directly regarded as being related to the ground level contribution. This processing yields a highly accurate forest height inversion for the RAMSES data [2].

Topographic fluctuations have been manifested to be one of the significant factors affecting the performance of the model-based Pol-InSAR inversion [6,27,28,29,30,31,32]. Specifically, the positive terrain slope causes an overestimation of the forest height, while the negative terrain slope causes an underestimation of the forest height [6,29,30,32]. For the purpose of eliminating the influence of terrain, the sloped Random Volume over Ground (S-RVoG) model was proposed by introducing the parameter of range terrain slope into the RVoG model [29]. With the addition of the model parameters, the multi-baseline optimization emerges as a direct inversion idea. Yet, when single-baseline optimization is still being implemented, prior terrain information is introduced to gain the terrain slope, which is then compensated into the RVoG model.

With the aim of obtaining a relatively precise volume coherence and simultaneously retaining the advantages of the single-baseline for inversion efficiency, this paper took advantage of the time-frequency optimization to separate ground and canopy contributions. The three-stage process is then modified with these two optimized coherences. Considering the influence of topographic fluctuations on forest height estimation, the S-RVoG model was adopted to improve the forest height inversion accuracy. The structure of this paper proceeds as follows. Section 2 presents the principle of the time-frequency optimization, while Section 3 introduces the S-RVoG model and illustrates the overall framework of model inversion. Afterwards, the verification of the real E-SAR data is presented in Section 4. Finally, Section 5 discusses and Section 6 concludes.

2. Time-Frequency Optimization Based on Sublook Decomposition

As mentioned earlier, to extract a relatively accurate volume coherence and to study the potential of the single-baseline configuration on P-band Pol-InSAR model inversion, the time-frequency optimization was proposed [2], which is based on synthetic aperture decomposition and Doppler filtering. This section will describe the principle and implementation of time-frequency optimization.

To improve the azimuth resolution, synthetic aperture radar (SAR) synthesizes a narrow beamwidth by time accumulation. The path that the sensor passes throughout this time is referred to as synthetic aperture length. During the synthetic aperture time, the Doppler history allows the azimuth signal to be treated as a chirp function, which is independent of the emitted signal. Therefore, as with signal processing in range, the azimuth focusing of the target imaging is also achieved by impulsion-compression. The top and center diagrams in Figure 1 indicate the synthetic aperture and Doppler history, respectively [42]. The synthetic aperture length pinpoints the actual antenna beamwidth, i.e., the scope of observation angle. In other words, the Doppler history is one-to-one with the radar observation angle. Therefore, the radar signal interrelating with the designated observation angle can be selected by filtering of Doppler spectra in the azimuth frequency-domain. Since diverse observation angles associate with diverse scattering characteristics [43,44,45], this time-frequency analysis can be utilized to separate the positions of various scattering centers [2,46].

The Doppler frequency [42] in Figure 1 can be expressed as a function of sensor speed v, radar wavelength $\lambda$, and observation angle ${\theta }_{sq}$

$$\begin{array}{c}\hfill f_D=\frac{2vsin{\theta }_{sq}}{\lambda }\end{array}$$

(1)

Since the Doppler frequency is proportional to the radial velocity of the target relative to the sensor, the absolute value of the Doppler frequency decreases first and then increases in the Doppler history. Furthermore, the Doppler frequency is positive when radar is approaching target, while presenting negative when the radar is moving away from the target. Therefore, the slope of the Doppler frequency versus time proves to be negative. The Doppler bandwidth [42] of the target is expressed as

$$\begin{array}{c}\hfill B_{Az}=\frac{2vcos{\theta }_{sqc}}{\lambda }\Delta {\psi }_{A{z}_{max}}\end{array}$$

(2)

where ${\theta }_{sqc}$ is the observation angle when the antenna beam center illuminates the target and $\Delta {\psi }_{A{z}_{max}}$ represents the 3-dB radar beamwidth. To fulfill the requirements of Nyquist sampling law, the bandwidth presented in Equation (2) determines the lower limit of the pulse-repetition frequency (PRF), i.e., $PRF\ge B_{Az}$. When PRF is greater than bandwidth, the normalized Doppler frequency spectra is specified in the bottom diagram of Figure 1 [42].

Marked with three different colors in Figure 1, the Doppler spectra is evenly divided into three parts through spectral filtering, overlapping each other by 50%. In accordance with Equation (1), the antenna patterns and observation angles corresponding to the three-segment spectra are also presented in Figure 1. Thus, the complete synthetic aperture is decomposed into three sublooks. Since a forest characterized by contrasted vegetation-density areas appears to have different scattering characteristics with various observation angles, the backscattered signal consistent with a specific observation angle allows to separate canopy and ground contributions, represented as high-density and low-density regions, respectively [2]. In the time-frequency analysis, the Doppler bandwidth scope is selected from the range-Doppler domain signal obtained by an azimuth Fourier transform, and then the sublook signal related to a specific observation angle is acquired employing an azimuth inverse Fourier transform. Since GVR varies with the observation angle for a given polarimetric channel [2], the separation of forest scattering centers is achieved by digging into the coherences of different sublooks. To sum up, the time-frequency optimization has the ability to independently extract backscattered signals at different observation angles, which means that it is possible to acquire a more accurate volume contribution. Therefore, this paper utilizes the optimization to improve the single-baseline P-band Pol-InSAR model inversion.

3. Model-Based Pol-InSAR Inversion

3.1. RVoG Model and Three-Stage Inversion Process

In the light of the assumption of the Volume over Ground [10,12] model, i.e., volume and ground contributions $S_v$ and $S_g$, the interferometric coherence between the master and slave antennas $S_m$ and $S_s$ for a given polarimetric state $\mathit{\omega }$ is derived as follows:

$$\begin{array}{c}\gamma \left(\mathbf{\omega }\right)=\frac{〈S_m\left(\mathbf{\omega }\right){S_s}^{\ast }\left(\mathbf{\omega }\right)〉}{\sqrt{〈S_m\left(\mathbf{\omega }\right){S_m}^{\ast }\left(\mathbf{\omega }\right)〉〈S_s\left(\mathbf{\omega }\right){S_s}^{\ast }\left(\mathbf{\omega }\right)〉}}\hfill \\ =\frac{〈\left(S_{g_m}\left(\mathbf{\omega }\right)+S_{v_m}\left(\mathbf{\omega }\right)\right){\left(S_{g_s}\left(\mathbf{\omega }\right)+S_{v_s}\left(\mathbf{\omega }\right)\right)}^{\ast }〉}{\sqrt{〈\left(S_{g_m}\left(\mathbf{\omega }\right)+S_{v_m}\left(\mathbf{\omega }\right)\right){\left(S_{g_m}\left(\mathbf{\omega }\right)+S_{v_m}\left(\mathbf{\omega }\right)\right)}^{\ast }〉〈\left(S_{g_s}\left(\mathbf{\omega }\right)+S_{v_s}\left(\mathbf{\omega }\right)\right){\left(S_{g_s}\left(\mathbf{\omega }\right)+S_{v_s}\left(\mathbf{\omega }\right)\right)}^{\ast }〉}}\hfill \\ =\frac{〈S_{g_m}\left(\mathbf{\omega }\right){S_{g_s}}^{\ast }\left(\mathbf{\omega }\right)〉+〈S_{v_m}\left(\mathbf{\omega }\right){S_{v_s}}^{\ast }\left(\mathbf{\omega }\right)〉}{{\left|S_g\left(\mathbf{\omega }\right)\right|}^2+{\left|S_v\left(\mathbf{\omega }\right)\right|}^2}\hfill \\ =\frac{\frac{〈S_{g_m}\left(\mathbf{\omega }\right){S_{g_s}}^{\ast }\left(\mathbf{\omega }\right)〉}{{\left|S_g\left(\mathbf{\omega }\right)\right|}^2}\frac{{\left|S_g\left(\mathbf{\omega }\right)\right|}^2}{{\left|S_v\left(\mathbf{\omega }\right)\right|}^2}+\frac{〈S_{v_m}\left(\mathbf{\omega }\right){S_{v_s}}^{\ast }\left(\mathbf{\omega }\right)〉}{{\left|S_v\left(\mathbf{\omega }\right)\right|}^2}}{\frac{{\left|S_g\left(\mathbf{\omega }\right)\right|}^2}{{\left|S_v\left(\mathbf{\omega }\right)\right|}^2}+1}\hfill \\ =e^{j{\varphi }_0}\frac{{\gamma }_g\mu \left(\mathbf{\omega }\right)+{\gamma }_v}{\mu \left(\mathbf{\omega }\right)+1}\hfill \end{array}$$

(3)

where the subscripts m and s represent the master and slave antennas, respectively; * indicates the conjugate operation; $〈〉$ means the expectation; and ${\varphi }_0$, ${\gamma }_g$, $\mu \left(\mathit{\omega }\right)$, and ${\gamma }_v$ denote the ground interferometric phase, ground coherence, GVR, and volume coherence, respectively. Notice that, in Equation (3), the ground and volume contributions are recognized as irrelevant. Interferometric coherence can be altered to the following form [10,12]

$$\begin{array}{c}\gamma \left(\mathbf{\omega }\right)=e^{j{\varphi }_0}\left[{\gamma }_v+L\left({\gamma }_g-{\gamma }_v\right)\right]\hfill \\ L\left(\mathbf{\omega }\right)=\frac{\mu \left(\mathbf{\omega }\right)}{1+\mu \left(\mathbf{\omega }\right)}\hfill \\ 0\le L\left(\mathbf{\omega }\right)<1\hfill \\ 0\le \mu \left(\mathbf{\omega }\right)<+\infty \hfill \end{array}$$

(4)

In general, ground coherence is assumed to be one, i.e., ${\gamma }_g=1$. It is concluded from Equation (4) that, for a given polarimetric state $\mathit{\omega }$, the ground coherence ${\gamma }_g$, interferometric coherence $e^{-j{\varphi }_0}\gamma \left(\mathit{\omega }\right)$, and volume coherence ${\gamma }_v$ are distributed on a straight line on the complex plane.

The RVoG [10,12] model adopting an exponential function to characterize the forest vertical backscatter profile deems that the volume coherence does not hinge on the polarimetric state, i.e., the assumption of random volume. On the basis of that assumption, the volume coherence with a forest height $h_v$ is represented as follows:

$$\begin{array}{c}\hfill {\gamma }_{RVoG}=\frac{{\int }_0^{h_v}{e}^{\frac{2\sigma z}{cos\theta }}{e}^{j{k}_{z}z}dz}{{\int }_0^{h_v}{e}^{\frac{2\sigma z}{cos\theta }}dz}\end{array}$$

(5)

where $\sigma$ and $\theta$ denote, respectively, the mean extinction coefficient and the radar incidence angle. The interferometric vertical wavenumber $k_z$ depends on the wavelength $\lambda$, incidence angle $\theta$, and its difference between the master and slave antennas $\Delta \theta$

$$\begin{array}{c}\hfill k_z=\frac{4\pi }{\lambda sin\theta }\Delta \theta \end{array}$$

(6)

As is shown in Figure 2, since ${\gamma }_{RVoG}$ is independent of the polarimetric state, when the GVR varies, all the homologous interferometric coherences calculated in light of Equation (4) move on a straight line on the complex plane. In view of that geometric characteristic, the single-baseline Pol-InSAR inversion is divided into three stages [12]:

(1)

Coherence line fit: According to the linear distribution prediction of the coherences in Equation (4) or Figure 2, the coherence line can be achieved by fitting the coherences of different polarimetric channels.

(2)

Ground interferometric phase estimation: Similarly, according to Equation (4) or Figure 2, the ground interferometric phase corresponds to one of the two intersections of the coherence line with the unit circle. The criteria of judging the true ground phase is as follows

$${\varphi }_0=\left\{\begin{array}{c}{\varphi }_{0_1},\left|e^{j{\varphi }_{0_1}}-{\gamma }_{HH+VV}\right|<\left|e^{j{\varphi }_{0_1}}-{\gamma }_{HV}\right|\hfill \\ {\varphi }_{0_2},else\hfill \end{array}\right.$$

(7)

where ${\varphi }_{0_1}$ and ${\varphi }_{0_2}$ indicate the two intersections of the coherence line with the unit circle. Equation (7) is based on the rank order of the polarimetric phase center positions, i.e., the HH+VV channel has a higher GVR than the HV channel, as shown in Figure 2. Since the phase diversity [47] can achieve a maximum phase separation, the PDHigh and PDLow channels can replace the HV and HH+VV channels in Equation (7).

(3)

Forest height estimation: Assuming that the GVR of a polarimetric channel (such as the PDHigh channel) is sufficiently low, in accordance with Equation (4), the volume coherence can be expressed as $e^{-j{\varphi }_0}{\gamma }_{PDHigh}$, and then the forest height is extracted by a two-dimensional look-up table (namely, the theoretical random volume solution space in Figure 2) set up in the light of Equation (5)

$$\begin{array}{c}\hfill {\widehat{h}}_v=\underset{\left(h_v,\sigma \right)}{Argmin}\left\{\left|e^{-j{\varphi }_0}{\gamma }_{PDHigh}-{\gamma }_{RVoG}\left(h_v,\sigma \right)\right|\right\}\end{array}$$

(8)

3.2. S-RVoG Model

Since the underlying topography over the forest area commonly takes on certain fluctuations, the forest height inversion on the basis of the traditional RVoG model under the flat ground hypothesis reckons to be susceptible to the terrain slope [6,27,28,29,30,31,32]. Accordingly, it is indispensable to incorporate the influence of the terrain slope into the forest height estimation. As shown in Figure 3, the terrain slope that tilted toward the radar called positive slope (see Figure 3a) decreases $\theta$, whereas the slope that tilted away from the radar called negative slope (see Figure 3b) increases $\theta$[30]. For reducing the complexity of inversion in consideration of the topography, the coordinate system is modified according to the terrain slope [29], denoted with $y{}^{\prime }o{}^{\prime }z{}^{\prime }$ in Figure 3. In the corrected coordinate system $y{}^{\prime }o{}^{\prime }z{}^{\prime }$, the three-stage method is continuously effective for forest height inversion. Neglecting the influence of azimuth, the terrain slope $\alpha$ is achieved by the following equation [30]

$$\begin{array}{c}\hfill \alpha ={tan}^{-1}\left(\frac{\Delta H}{\frac{P{S}_{rg}}{sin\theta }+\frac{\Delta H}{tan\theta }}\right)\end{array}$$

(9)

where $\Delta H$ calculated from the measured DEM is the altitude difference between two consecutive slant range samples and ${PS}_{rg}$ represents the homologous slant range distance. In Figure 3, the terrain slope $\alpha$ corresponding to a positive slope takes a positive value, while that corresponding to a negative slope takes a negative value. Thus, the local incidence angle corrected by $\alpha$ in the corrected coordinate system $y{}^{\prime }o{}^{\prime }z{}^{\prime }$ is then rewritten as

$$\begin{array}{c}\hfill \theta {}^{\prime }=\theta -\alpha \end{array}$$

(10)

Since the local incidence angle is corrected by the terrain slope, interferometric vertical wavenumber $k_z$ also demands to be adjusted accordingly, while $\Delta \theta$ remains unchanged before and after the coordinate system correction. Thus, the corrected vertical wavenumber can be expressed as

$$\begin{array}{c}\hfill k_z{}^{\prime }=\frac{4\pi }{\lambda sin\theta {}^{\prime }}\Delta \theta \end{array}$$

(11)

Equations (5) and (6) present that $\theta$ and $k_z$ in volume coherence are directly dependent on the change in incidence angle. Therefore, in the corrected coordinate system $y{}^{\prime }o{}^{\prime }z{}^{\prime }$, the homologous volume coherence on account of the S-RVoG model requires an adjustment according as Equations (10) and (11), as shown in the following formula

$$\begin{array}{c}\hfill {\gamma }_{S-RVoG}=\frac{{\int }_0^{h_v{}^{\prime }}{e}^{\frac{2\sigma z{}^{\prime }}{cos\theta {}^{\prime }}}{e}^{j{k}_z{}^{\prime }z{}^{\prime }}dz{}^{\prime }}{{\int }_0^{h_v{}^{\prime }}{e}^{\frac{2\sigma z{}^{\prime }}{cos\theta {}^{\prime }}}dz{}^{\prime }}\end{array}$$

(12)

where $h_v{}^{\prime }$ indicates the forest height in coordinate system $y{}^{\prime }o{}^{\prime }z{}^{\prime }$. It is concluded from Equations (10)–(12) that the positive slope increases $k_z$, resulting in a decrease in the volume coherence, while the negative slope produces an opposite effect [30,32]. The essence of the volume coherence represented by Equation (12) is a process in which the accumulation of phase changes caused by $k_z{}^{\prime }$ is normalized. It manifests that when $k_z{}^{\prime }$ is equal to zero, namely the phase remains at zero, the volume coherence calculated in accordance with Equation (12) is identically equal to one. Under the circumstance, the calculation of volume coherence is independent of $\theta {}^{\prime }$. The volume coherence reaches the maximum value at this time. When $k_z{}^{\prime }$ is greater than zero, the local incidence angle $\theta {}^{\prime }$ will also have an impact on the volume coherence. However, from Equation (12), $\theta {}^{\prime }$ has less effect than that of $k_z{}^{\prime }$ on the volume coherence. Consequently, as $k_z{}^{\prime }$ becomes larger, the volume coherence becomes smaller [30].

Since $y{}^{\prime }o{}^{\prime }z{}^{\prime }$ is projected from $yoz$, the forest height estimation $h_v{}^{\prime }$ and $h_v$ in the two coordinate systems equally have a projection relationship, as shown in Figure 3. Hence, $h_v$ perpendicular to the horizon is obtained by performing an inverse projection operation on $h_v{}^{\prime }$. The pertinent expression is as follows

$$\begin{array}{c}\hfill h_v=\frac{h_v{}^{\prime }}{cos\left|\alpha \right|}\end{array}$$

(13)

In addition, as can be seen from Figure 3, the relative extinction path has also changed in the corrected coordinate system, which directly affects the GVR. Therefore, Equation (4) is adjusted to the following form

$$\begin{array}{c}\gamma {}^{\prime }\left(\mathbf{\omega }\right)=e^{j{\varphi }_0}\left[{\gamma }_v{}^{\prime }+L{}^{\prime }\left({\gamma }_g-{\gamma }_v{}^{\prime }\right)\right]\hfill \\ L{}^{\prime }\left(\mathbf{\omega }\right)=\frac{\mu {}^{\prime }\left(\mathbf{\omega }\right)}{1+\mu {}^{\prime }\left(\mathbf{\omega }\right)}\hfill \\ 0\le L{}^{\prime }\left(\mathbf{\omega }\right)<1\hfill \\ 0\le \mu {}^{\prime }\left(\mathbf{\omega }\right)<+\infty \hfill \end{array}$$

(14)

where the superscript ${}^{\prime }$ represents the parameters in the corrected coordinate system.

The concrete effect of the terrain slope on the volume coherence is shown in Figure 4. Compared with the random volume solution space in Figure 2, the volume coherence decreases on a positive terrain slope (see Figure 4a) while increases on a negative terrain slope (see Figure 4b), manifesting a consistence with the previous theoretical analysis. The pink and green curves represent ${\gamma }_{S-RVoG}$ and ${\gamma }_{RVoG}$ for the same forest height, respectively. The findings reveal that, with the random volume solution space inversion, the positive terrain slope will cause an overestimation (see Figure 4a), whereas the negative terrain slope will cause an underestimation (see Figure 4b). Figure 4 presents the result after the projection relationship expressed by Equation (13) has been merged into the volume coherence calculated in line with Equation (12).

3.3. S-RVoG Model Inversion Based on Time-Frequency Optimization

Techniques such as polarimetric synthesis and coherence optimization [10,47] select different scattering centers from the scattering mechanism dimension, while time-frequency optimization separates phase centers from the observation angle dimension. For the purpose of achieving a maximum separation of the vertical phase center, in this paper, the time-frequency optimization is adopted to perform sublook decomposition imaging and then the highest (TFHigh) and lowest (TFLow) phase centers can be derived from the complex coherences of the selected polarimetric channels, which are estimated separately for each sublook image pair. To mitigate the effects of temporal decorrelation and noise, phase diversity [47] and magnitude diversity [10] optimizations are employed to pick the polarimetric channels for coherences estimation [48]. In regard to three-stage method, the polarimetric channel with a sufficiently small GVR [12] is pinpointed to estimate the volume coherence. Since the combination of the time-frequency analysis and coherence optimization can fulfill an effective separation of the scattering centers, it is reasonable to assume that a TFHigh channel corresponds to a relatively small GVR.

The S-RVoG inversion process is shown in Figure 5, in which ${\gamma }_{TFH\mathrm{igh}}$ is presumed to be close enough to the real volume coherence [2]. Possibly, there is still some distance between ${\gamma }_{TFH\mathrm{igh}}$ and the real volume coherence. The green curve indicates the theoretical sloped random volume for 0 dB/m extinction, namely the boundary of the solution space. When ${\gamma }_{TFH\mathrm{igh}}$ does not fall into the solution space, the minimum Euclidean distance criterion referenced in the look-up table method provides a predicted volume coherence on the zero extinction curve which is represented in the blue dot in Figure 5. After preprocessing the data with the time-frequency optimization, the inversion process of the S-RVoG model is as follows. It should be noted that the process is still rooted on the framework of the traditional three-stage method.

(1)

Coherence line fit: This stage is consistent with the traditional three-stage approach. The total least squares method is adopted to fit the coherence line. ${\gamma }_{TFH\mathrm{igh}}$ and ${\gamma }_{TFL\mathrm{ow}}$ will take part in the fitting of the coherence line, further suppressing the effects of temporal decorrelation and noise.

(2)

Ground interferometric phase estimation: The ground interferometric phase is still acquired by one of the intersections of the coherence line and the unit circle, yet the criteria for selecting the real ground interferometric phase is adjusted

$${\varphi }_0=\left\{\begin{array}{c}{\varphi }_{0_1},\left|e^{j{\varphi }_{0_1}}-{\gamma }_{TFLow}\right|<\left|e^{j{\varphi }_{0_1}}-{\gamma }_{TFHigh}\right|\hfill \\ {\varphi }_{0_2},else\hfill \end{array}\right.$$

(15)

(3)

Forest height estimation: Since the GVR related to ${\gamma }_{TFH\mathrm{igh}}$ is reckoned to be sufficiently small, that is, homologous $L{}^{\prime }$ in Equation (14) is approximately equal to zero, the volume coherence is expressed as $e{}^{-j{\varphi }_0}{\gamma }_{TFHigh}$. Therefore, the forest height is estimated from a two-dimensional look-up table which is generated in the light of Equation (12) on the basis of the minimum Euclidean distance criterion

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\widehat{h}}_v{}^{\prime }=\underset{\left(h_v{}^{\prime },\sigma \right)}{Argmin}\left\{\left|e^{-j{\varphi }_0}{\gamma }_{TFHigh}-{\gamma }_{S-RVoG}\left(h_v{}^{\prime },\sigma \right)\right|\right\}\\ {\widehat{h}}_v=\frac{{\widehat{h}}_v{}^{\prime }}{cos\left|\alpha \right|}\end{array}\right.\end{array}$$

(16)

where ${\widehat{h}}_v{}^{\prime }$ denotes the estimator of $h_v{}^{\prime }$. Figure 6 portrays the flowchart of forest height inversion based on the time-frequency optimization reckoning terrain compensation. The blue block diagrams represent the key processing steps in this paper.

4. Forest Height Estimation Using P-Band E-SAR Data

4.1. Test Site and E-SAR Data Description

This paper performs a method validation relying on the European Space Agency (ESA) BIOSAR 2008 data [6] acquired by the German Aerospace Center (DLR) E-SAR system in Krycklan catchment, Northern Sweden. The test site is mainly covered by boreal forest that is primarily between 0–30 m in height. The data set contains the lidar-derived forest height, which is applied to evaluate the performance of the inversion method. In addition, the measured topography elevation of the experimental area is utilized in this paper to generate terrain slope and to compensate for the inversion errors caused by topographic fluctuations. Figure 7 presents a pauli–basis composite map, measured DEM, and lidar-derived forest height of the test site.

Table 1. The parameters of the interferometric pair.

Track	Baseline (m)	${\mathit{k}}_{\mathit{z}}$ Range	Band	Polarization
Master	Master	Master	P	Quad
Slave	32	0.026–0.245	P	Quad

shows the single-baseline data information researched in this paper.

4.2. Forest Height Estimation

In accordance with the flowchart represented in Figure 6, the sub-aperture imaging under four sublook decompositions is first performed along the azimuth direction, and then, the two sets of sublook images are subjected to interferometric processing to estimate the coherences of the specified polarimetric channels. The phase diversity [47] applied in this process further optimizes the coherences from the perspective of scattering mechanism. ${\gamma }_{TFH\mathrm{igh}}$ and ${\gamma }_{TFL\mathrm{ow}}$ are employed in the modified three-stage algorithm to improve the failure of the algorithm at P-band observations. Simultaneously, the measured DEM providing the range terrain slope is utilized to correct the local incidence angle and interferometric vertical wavenumber (see Equations (10) and (11)). Here, the effects of azimuth terrain fluctuations are not taken into account. In the corrected coordinate system, the sloped random volume solution space calculated in the light of the theoretical volume coherence is also updated (see Equation (12)), and based on that, a two-dimensional look-up table method is adopted to extract the forest height.

In order to separately test the inversion performance improved by time-frequency optimization and terrain compensation, four sets of comparative experiments are each implemented. Concretely, four experimental configurations are carried out, i.e., (a) RVoG model + three-stage method, (b) RVoG model + modified three-stage method (time-frequency optimization), (c) S-RVoG model + three-stage method, and (d) S-RVoG model + modified three-stage method (time-frequency optimization). The corresponding inversion maps are presented in Figure 8, illustrating that the inversions of the four cases manifest an identical spatial trend and that the traditional three-stage method generates a serious underestimation, while the time-frequency optimization can effectively ameliorate this underestimation. Note that, although continuing to increase the number of sublook will further improve the inversion performance, the reduction in the azimuth resolution is another limit worth considering. By comparing the difference value maps after terrain compensation (see Figure 8e,f) and the slope map (see Figure 8g), a conclusion consistent with the theoretical analysis is obtained, that is, the positive slope compensation reduces the estimation while the negative slope compensation raises the estimation. Further, quantify the four inversions in Figure 9. With the aim of accurately evaluating the performance, 200 validation stands are randomly selected to quantify the forest height inversion accuracy, as shown in Figure 9. It can be seen from Figure 8 or Figure 9 that it is necessary to perform the time-frequency optimization processing before terrain compensation. As a result, in the following analysis, this paper focuses on three cases, as shown in Figure 8a,b,d or Figure 9a,b,d. It is reflected from Figure 9a,b,d that the time-frequency optimization improves the inversion performance by 41.49%, while the improvement by terrain compensation is relatively insignificant on the whole. From the previous analysis, the influence of the positive and negative terrain slopes on the inversion performance is not necessarily symmetrical (see Figure 4), thus a comparative analysis is presented in the two terrain slope regions respectively, as shown in Figure 10 (corresponding to Figure 9b,d). The results manifest that the terrain compensation on the positive terrain slope decreases the inversion performance by 8.43%. According to the theoretical analysis, when inversion is on the basis of the RVoG model, the positive terrain slope will generally cause an overestimation, yet for the P-band data, this configuration presents an underestimation on the whole (see Figure 10a) due to the inaccuracy of the volume coherence estimation. Since the S-RVoG model on the positive terrain slope region will lower the inversion (see Figure 4a), the region will receive a more serious underestimation after the terrain compensation, which is the reason for the performance degradation of the region. With regard to the negative terrain slope region, the root mean square error (RMSE) depending on the RVoG model is 5.35 m, whereas that depending on the S-RVoG model is 4.71 m, which improves the performance by 11.96%, manifesting a consistence with the theoretical analysis. For the sake of presenting more intuitively the effect of terrain compensation on positive and negative slope regions, density maps of the inversion error with slope distribution based on the RVoG and S-RVoG models are shown in Figure 10g,h, respectively. It is significant that terrain compensation reduces the inversion on the positive terrain slope, while raises the inversion on the negative terrain slope. In fact, if there is a regular deviation between inversion and measurement, the deviation can be easily corrected by ancillary data. However, as can be seen from Figure 10i, the data researched in this paper is not suitable for that correction strategy.

5. Discussion

5.1. Phase Center Separation Based on Time-Frequency Optimization

The essence of synthetic aperture is the accumulation of the azimuth time, enabling us to study the radar echo signals related to specific observation angle gaps by selecting the Doppler frequency. Accordingly, the time-frequency analysis enables us to separate the interferometric phase centers from the dimension of the radar observation angle, as different observation angles generally illuminate various scattering mechanisms. For high-frequency SAR systems, the HV channel is commonly deemed to be close to the channel associated with volume scattering, noting that the HV channel does not necessarily correspond directly to a pure volume contribution. Nevertheless, this assumption is no longer valid under P-band observations. As is presented in Figure 11, the HV channel has been quite different from the pure volume channel. Although the coherence optimization acquires a channel, namely ${\gamma }_{PDH\mathrm{igh}}$, of which the phase center is higher than that of the HV channel, the time-frequency optimization further effectively separates the phase centers, achieving a higher one (see ${\gamma }_{TFH\mathrm{igh}}$ in Figure 11). Simultaneously, the phase center of ${\gamma }_{TFL\mathrm{ow}}$ gained with the time-frequency optimization is also significantly lower than that of ${\gamma }_{PDL\mathrm{ow}}$ gained with the coherence optimization. Accordingly, the effectiveness of a time-frequency optimization can be directly verified by Figure 11. The result shown in Figure 11 corresponds to the four sublooks situation. As the number of sublook increases, the phase center will be further separated; nonetheless, the loss of the azimuth resolution will become greater [2].

5.2. Effects of Terrain Slope on Inversion Sensitivity

In accordance with Equations (10) and (11), the positive slope increases $k_z$ while the negative slope decreases $k_z$. In the light of Equation (12), a larger $k_z$ leads to the decrease of coherence and the increase of the phase in the random volume solution space, while a smaller $k_z$ leads to an opposite effect. That conclusion is also verified by comparing Figure 2 and Figure 4, from which we can also conclude that the solution space corresponding to a larger $k_z$ owns a lower inversion sensitivity, whereas that corresponding to a smaller $k_z$ owns a higher inversion sensitivity [6,35]. The homologous experimental results are presented in Figure 12. Consistent with the theoretical analysis, the S-RVoG model on a positive slope reduces the inversion, whereas that on a negative slope rises. Compared with the RVoG model inversion, the S-RVoG model with a negative slope possesses a significantly greater inversion change than that with positive slope, which confirms the conclusion that the negative slope owns a greater inversion sensitivity. In Figure 12a, there are several forest stands of which the inversion changes are located in the scope larger than zero. That accounts for a larger $k_z$ comprising a small height of ambiguity. To remove these invalid pixels, we have interpolated these pixels, which brings about the occurrence of anomalous stands in Figure 12a.

5.3. Discussion on the Quality of Forest Height Estimation

Various forest height estimation methods have been proposed for the P-band data of BIOSAR 2008 campaign [6,32,35,48]. Whether it uses Pol-InSAR or TomoSAR to retrieve forest height, the multi-baseline configuration is a very important solution. In this case, the scheme still using single-baseline configuration requires some prior information. Compared with these studies, the time-frequency optimization is achieved under the single-baseline configuration and does not require specific prior information. Correspondingly, while improving the inversion efficiency, it loses some inversion accuracy. Even so, from the published literature, the performance of this paper is very close to the results of some multi-baseline configurations (such as dual-baseline). To sum up, the proposed method is applicable to the situation that the multi-baseline acquisition cannot be satisfied, there is no specific prior information, a high inversion efficiency is required, and a certain loss of inversion accuracy can be tolerated.

5.4. Limitations of the Proposed Method

From the results discussed in Section 4, even if the time-frequency optimization ameliorates the underestimation of the traditional three-stage method, the underestimation phenomenon still exists. That accounts for at P-band observations all polarimetric channels involving a certain ground contribution. Future research can try to combine the multi-baseline configuration with the time-frequency optimization to avoid the NGVR hypothesis. Another limitation is that the sublook decomposition loses the azimuth resolution of SAR image, and its quantification requires carrying out in future research. The last limitation is that the RVoG model may not be suitable for depicting the vertical backscatter profile of P-band Pol-InSAR as the time-frequency method still cannot completely eliminate the underestimation. That is because the predicted volume coherence is at the boundary (i.e., zero attenuation) of the solution space, which indicates that the situation is beyond the scope of the RVoG model description. The RVoG model deems that the power peak locates at the top of the canopy, nevertheless, when P-band is utilized to observe the boreal forest, backscatter power is chiefly concentrated in the middle and lower parts of the canopy. The way to cope with this limitation is to adopt some new complicated models [35,48,49,50,51].

6. Conclusions

In this paper, we choose the single-baseline configuration for the Pol-InSAR model inversion. Since P-band observations bring about a serious underestimation of the traditional three-stage method, the time-frequency optimization is taken advantage of to mitigate the phenomenon. The phase centers are effectively separated from the dimension of the observation angle, which increases the inversion performance by 41.49%. Furthermore, we take into account the influence of the range terrain slope on inversion, compensating for the three-stage method based on the time-frequency optimization, which improves the inversion performance on negative terrain slope by 11.96%. In summary, this paper provides a scheme on the basis of a time-frequency analysis to solve the Pol-InSAR inversion problem of the single-baseline P-band configuration for sparse vegetation observations. Simultaneously, terrain correction processing is effectively incorporated into this scheme.