Forest Height Inversion Based on Time–Frequency RVoG Model Using Single-Baseline L-Band Sublook-InSAR Data

Abstract

The interferometric synthetic aperture radar (InSAR) technique based on time–frequency (TF) analysis has great potential for mapping the forest canopy height model (CHM) at regional and global scales, as it benefits from the additional InSAR observations provided by the sublook decomposition. Meanwhile, due to the wider swath and higher spatial resolution of single-polarization data, InSAR has a higher observation efficiency in comparison with PolInSAR. However, the accuracy of the CHM inversion obtained by the TF-InSAR method is attenuated by its inaccurate coherent scattering modeling and uncertain parameter calculation. Hence, a new approach for CHM estimation based on single-baseline InSAR data and sublook decomposition is proposed in this study. With its derivation of the coherent scattering modeling based on the scattering matrix of sublook observations, a time–frequency based random volume over ground (TF-RVoG) model is proposed to describe the relationship between the sublook coherence and the forest biophysical parameters. Then, a modified three-stage method based on the TF-RVoG model is used for CHM retrieval. Finally, the two-dimensional (2-D) ambiguous error of pure volume coherence caused by residual ground scattering and temporal decorrelation is alleviated in the complex unit circle. The performance of the proposed method was tested with airborne L-band E-SAR data at the Krycklan test site in Northern Sweden. Results show that the modified three-stage method provides a root-mean-square error (RMSE) of 5.61 m using InSAR and 14.3% improvement over the PolInSAR technique with respect to the classical three-stage inversion result. An inversion accuracy of RMSE = 2.54 m is obtained when the spatial heterogeneity of CHM is considered using the proposed method, demonstrating a noticeable improvement of 32.8% compared with results from the existing method which introduces the fixed temporal decorrelation factor.

1. Introduction

A forest canopy height model (CHM) is important for the detection of forest harvesting, forest degradation, and forest fires [1]. Moreover, timely access to accurate estimates of the above-ground forest biomass would assist in studying the regional and global carbon cycle and climatic variation [2,3,4,5].

At present, ground measurement, light detection and ranging (LiDAR), and interferometric synthetic aperture radar (InSAR) have been widely applied to estimate the CHM for different forest types (boreal, temperate, and tropical forests). Due to the low efficiency and the high cost of data acquisition, it is difficult for ground measurement to meet the requirements of rapid and large-scale CHM mapping. Meanwhile, the accuracy of the measured forest CHM is limited by the density of sampling points [6]. LiDAR and InSAR, as active remote sensing techniques, have significant advantages over ground-based measurements in terms of the scope of data observation. However, the application of ground-based or airborne LiDAR is limited by the coverage of data and the exorbitant price [7,8]. In addition, the spaceborne LiDAR (e.g., ICESat [9]) is susceptible to cloud interference and has difficulty penetrating dense forests. Fortunately, InSAR has great potential for large-scale and high-accuracy CHM mapping because of the strong penetration ability of its radar wave, its independence from weather conditions, and its sensitivity to the vertical distribution of scatterers [10].

The CHM inversion methods based on the InSAR technique can be divided into three categories, which are described as follows.

(1) Traditional InSAR inversion: There are two ways to estimate the CHM from InSAR observations. On the one hand, the CHM can be extracted from the gap between the center height of the interference phase of the single-polarization InSAR and the external underlying DEM [11,12]. It is inevitable that the accuracy of the CHM estimates largely depends on the high-precision of the external underlying DEM, which is generally difficult to obtain. On the other hand, the degraded random volume over ground (RVoG) model can also be used to estimate the CHM from InSAR data, but it reduces the accuracy of the CHM.

(2) PolInSAR inversion: Polarimetric SAR interferometry (PolInSAR) is developed from InSAR by combining it with multi-polarization information. With this technique, different scatterers which are located in the vertical dimension of the same pixel can be detected, and then the CHM estimates can be obtained without using other external data [13,14,15,16]. Moreover, the multi-baseline technique developed in recent years further provides a variety of spatial observation schemes for PolInSAR to improve the accuracy of the CHM inversion [17,18,19,20]. As a result, these improved techniques can be used to avoid the underdetermined problem of unknown parameters based on a vegetation scattering model (e.g., RVoG), which makes it possible to obtain the high-precision CHM. However, the CHM inversion based on the multi-polarization data comes at the cost of reducing the observation efficiency of InSAR. This is because the PolInSAR technique usually has a smaller swath and lower spatial resolution than the single-polarization case [21].

(3) Sublook InSAR inversion: The time-frequency (TF) analysis method provides a new pathway for high-precision CHM inversion. In previous studies, the dual-polarization TF optimization method has been used to estimate the accurate interference phase of pure volume coherence and the ground phase [22], which helps to improve the accuracy of the CHM inversion. Besides, the TF optimization method has also been extended to the full-polarization case [23]. It can be seen that the TF analysis can be used as an improved method to obtain the optimal coherence. However, compared with single-polarization InSAR, the dual- and full-polarization cases will lower the observation efficiency.

In the single-polarization case, the Time-frequency and Line-fitting (TF + LF) method with a priori extinction coefficient has been proposed for the CHM inversion from the single-baseline InSAR data. It can not only solve the underdetermined problem of unknown parameters but also avoid the low observation efficiency of the PolInSAR technique [6,24]. However, the performance of this method is usually affected by three factors. The first one comes from the coherent scattering model. In the TF + LF method, the sublook coherence plays a similar role to that of the RVoG-based polarization coherence. However, an essential difference between the sublook and polarimetric observations has not been revealed in the TF + LF method, which leads to an ambiguous interpretation of the impact of the sublook observations. As for the difference, it is known that the polarization coherences are determined by the different polarization states, whereas the sublook coherences are distinguished by the path differences of the sublook SAR signals, which are caused by the varying azimuth observation angles. The second problem is referred to as the defective parameter inversion. Namely, a small amount of CHM is necessary to extract the extinction coefficient which will be used as a priori knowledge. The third problem is the impact of the pure volume coherence error. The estimated CHM is sensitive to the 2-D ambiguous error of pure volume coherence, which is induced by the joint impacts of the residual ground scattering and the temporal decorrelation and is difficult to removed in the single-baseline inversion.

Therefore, the CHM inversion method based on TF analysis has great potential to extract a high-precision CHM benefitting from the more available observations provided by the sublook decomposition and the higher observation efficiency of the single-polarization InSAR. However, the inaccurate coherent scattering modeling and parameter inversion are major problems that limit the accuracy of the CHM estimation.

According to the above methods, the information from each scatterer is collected at various incidence angles, and the sublook coherences along the azimuth carry different backscattering properties. Note that the essential reason for this is the path difference of the sublook SAR signals through the forest volume, which is important for the model-based CHM inversion from the sublook coherences, and should be further studied. As a result, the scattering matrix was reconstructed in consideration of the path differences of sublook observations in this paper. With a derivation of coherent scattering modeling based on the reconstructed scattering matrix, a more accurate sublook coherent scattering model, namely, the TF-RVoG model, is proposed for the CHM inversion from sublook InSAR data. Regarding the height of the phase center for different sublook coherences, the pure volume coherence can be obtained from the alternative sublook coherences. Hence, a single-baseline three-stage method for the TF-RVoG model was also proposed in this work. In order to obtain a higher accuracy for the CHM, the empirical relationship between $D.I$ (which is the distance ratio index used to describe the position of pure volume coherence in the coherence line [25]) and the extinction coefficient is used to alleviate the residual ground scattering error of the pure volume coherence. In addition, the strong correlation between the temporal decorrelation and the CHM is verified, and a random volume over ground with mutative temporal decorrelation (RVoG + MTD) model that takes into account the spatial heterogeneity of the forest CHM is proposed to remove the temporal decorrelation error. The advantages of the proposed approach are as follows:

(1)

It considered the path difference of sublook SAR signals and can help to provide detailed interpretation for the impact of the different sublook coherences in the sublook coherent scattering modeling;

(2)

It can invert the forest CHM from the single-baseline and single-polarization InSAR data without an a priori CHM, underdetermined parameter inversion, and observation efficiency reduction;

(3)

It can alleviate the influence of the 2-D ambiguous error of pure volume coherence by combining the empirical relationship of the $D.I$ and the extinction coefficient [25] with the RVoG + MTD model.

The rest of this paper is structured as follows. Section 2 introduces the principle and modeling process of the proposed TF-RVoG model, the correction of interferometric phase error of sublook coherences induced by the residual motion error (RME) [26], the three-stage CHM inversion process based on the TF-RVoG model, and the correction of 2-D ambiguous error of pure volume coherence. In Section 3, the performance of the proposed approach is tested using single-baseline L-band InSAR data acquired by the German Aerospace Center (DLR)’s airborne E-SAR system. Discussion of and conclusions about the proposed approach are presented in Section 4 and Section 5, respectively.

2. Theory and Method

2.1. Sublook Coherent Scattering Modeling

In the TF analysis, fast Fourier transform (FFT) filtering is used to split the full-resolution SAR image into several sublook images [27,28,29], and then the target can be detected using various azimuth observation angles from the different sublook SAR images. The azimuth observation angle can be expressed as a function of the Doppler frequency $f_D$, the radar wavelength $\lambda$, and the platform velocity $v_p$ [30]

$${\alpha }_{Az}=\arcsin \left(\frac{f_D\lambda }{2v_p}\right)$$

(1)

It is the angle between the sight line associated with a zero Doppler centroid and the random sight line, which determines the relative position of a sublook along the azimuth. Moreover,${\alpha }_{Az}$ can range from ${\alpha }_{Az}-0.5\times \Delta {\alpha }_{Az}$ to ${\alpha }_{Az}+0.5\times \Delta {\alpha }_{Az}$, and $\Delta {\alpha }_{Az}$ is the variation interval of the azimuth observation angle related to the azimuth resolution of the sublook image. As shown in Figure 1 (top), FFT filtering of the full-resolution SAR-image frequency spectra is performed in the azimuth direction to form three sublooks of equivalent resolution with a 50% overlap, which are represented as green, red, and blue dashed rectangles. Note that only the Doppler spectra of the red and blue sublooks are independent. Here, the green, red, and blue arrows in Figure 1 (bottom) denoted the azimuth observation angles which are related to the three sublooks. Since the filtered bandwidths are constant, the intervals of the observation angle are equivalent for the three sublooks.

Two sublook images $s_1$ and $s_2$ with the same ${\alpha }_{Az}$ and $\Delta {\alpha }_{Az}$ can be acquired by performing the sublook decomposition of the full-resolution master and slave images, respectively. After interferometric processing, the corresponding sublook coherence ${\tilde{\gamma }}_{Obs}$ can be expressed as [24]

$${\tilde{\gamma }}_{Obs}\left({\alpha }_{Az},\Delta {\alpha }_{Az}\right)=\frac{\langle s_1\left({\alpha }_{Az},\Delta {\alpha }_{Az}\right)\times s_2^{*}\left({\alpha }_{Az},\Delta {\alpha }_{Az}\right)\rangle }{\sqrt{\langle s_1\left({\alpha }_{Az},\Delta {\alpha }_{Az}\right)\times s_1^{*}\left({\alpha }_{Az},\Delta {\alpha }_{Az}\right)\rangle \langle s_2\left({\alpha }_{Az},\Delta {\alpha }_{Az}\right)\times s_2^{*}\left({\alpha }_{Az},\Delta {\alpha }_{Az}\right)\rangle }}$$

(2)

where 〈 〉 represents the expected value. ${\tilde{\gamma }}_{Obs}$ is impacted by several decorrelation factors, such as the signal-to-noise decorrelation, the geometric decorrelation, and the temporal decorrelation, which can be reduced or absorbed by the associated correction methods [31,32,33]. On this basis, the sublook coherence ${\tilde{\gamma }}_{Obs}$ is approximately equivalent to the volume decorrelation related to the backscattering properties of scatterers in forest scenes.

It is clear that the sublook coherences help to increase InSAR observations, which makes it possible to invert the CHM in the case of single-polarization InSAR. Therefore, a sublook coherent scattering model, which can be used to accurately describe the relationship between the sublook coherence and the forest biophysical parameters, must be built for the CHM inversion. However, the varying azimuth observation angles of the sublook coherences result in the path difference of the sublook signals in the forest volume, which can cause backscattering variations because of the anisotropic geometrical structures of the real forest volume. Therefore, the path impacts of the different sublook signals need to be considered and accounted for in the coherent scattering model.

Based on the above analysis, the scattering matrix can be expressed as [32]

$$S_{Obs}=WS{W}^T=\left[\begin{array}{cc}{e}^{-i{\beta }_{h}z}& 0\\ 0& e^{-i{\beta }_{v}z}\end{array}\right]·\left[\begin{array}{cc}{S}_{hh}& S_{hv}\\ S_{vh}& S_{vv}\end{array}\right]·\left[\begin{array}{cc}{e}^{-i{\beta }_{h}z}& 0\\ 0& e^{-i{\beta }_{v}z}\end{array}\right]$$

(3)

where the superscript $T$ denotes the transpose of a given matrix, and the impacts of the different propagation paths of sublook signals can be expressed with a 2 × 2 propagation matrix $W$ which is defined as

$$W=\left[\begin{array}{cc}{e}^{-i{\beta }_{h}z}& 0\\ 0& e^{-i{\beta }_{v}z}\end{array}\right]$$

(4)

where ${\beta }_h$ and ${\beta }_v$ are in general complex, relating to both phase shift and attenuation in the forest volume. Specifically, ${\beta }_h={\beta }_0{n}_h$, ${\beta }_v={\beta }_0{n}_v$, and ${\beta }_0$ is the free space wavenumber, which is equal to $2\pi /\lambda$. $n_h$ and $n_v$ represent the propagation value of the horizontal polarization and vertical polarization, respectively. The expression $S_{Obs}=WS{W}^T$ describes the process of the sublook signals’ propagation into and then out of the medium.

After vectorization processing based on the Pauli basis, the scattering matrix $S_{Obs}$ can be transformed into a four-dimensional (4-D) vector expressed as

$$\begin{gathered} K_{P_4}\left(z\right)=P_{P_4}·k_{P_4} \\ {P}_{P_4}=e^{-i\left({\beta }_h+{\beta }_v\right)z}\left[\begin{array}{cccc}\cos \mathrm{h}\left(sz\right)& \sin \mathrm{h}\left(sz\right)& 0& 0\\ \sin \mathrm{h}\left(sz\right)& \cos \mathrm{h}\left(sz\right)& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right] \\ \mathrm{where}\{\begin{array}{c}s=-i\left({\beta }_h-{\beta }_v\right);\\ \sinh \left(sz\right)=\frac{e^{sz}-e^{-sz}}{2};\\ \cosh \left(sz\right)=\frac{e^{sz}+e^{-sz}}{2};\end{array} \\ {k}_{p_4}=\frac{1}{2}{\left[\begin{array}{cc}\begin{array}{cc}{S}_{hh}+S_{vv}& S_{hh}-S_{vv}\end{array}& \begin{array}{cc}{S}_{hv}+S_{vh}& i\left(S_{hv}-S_{vh}\right)\end{array}\end{array}\right]}^T \end{gathered}$$

(5)

where $s$ is the differential propagation constant; $P_{P_4}$, as a 4 × 4 matrix, represents the propagation matrix associated with propagation path within the medium; and sinh and cosh are hyperbolic sine and hyperbolic cosine functions, respectively. In the case of a random medium, the parameter s is equal to zero, and $P_{P_4}$ reduces to a multiple of the identity matrix. Furthermore, $k_{P_4}$ is a 4-D scattering vector that does not the propagation influence on the vegetation backscattering. Based on the reciprocity theorem for antennae, namely$S_{hv}=S_{vh}$, the scattering vector $k_{P_4}$ can be simplified as a three-dimensional (3-D) vector expressed as

$$\begin{gathered} K_{P_3}\left(z\right)=P_{P_3}·k_{P_3} \\ {P}_{P_3}=e^{-i\left({\beta }_h+{\beta }_v\right)z}\left[\begin{array}{ccc}cosh\left(sz\right)& sinh\left(sz\right)& 0\\ sinh\left(sz\right)& cosh\left(sz\right)& 0\\ 0& 0& 1\end{array}\right] \\ {k}_{P_3}=\frac{1}{2}{\left[\begin{array}{cc}\begin{array}{cc}{S}_{hh}+S_{vv}& S_{hh}-S_{vv}\end{array}& 2S_{hv}\end{array}\right]}^T \end{gathered}$$

(6)

where $P_{P_3}$ is the corresponding 3 × 3 propagation matrix, and $k_{P_3}$ is the scattering vector which is free from the propagation influence.

According to the three-component decomposition based on the different scattering mechanisms, the scattering vector $P_{P_3}$ can be expressed as

$$\begin{gathered} K_{P_1}=P_{P_3}·k_{P_3-master}=k_{v1}+P_s{k}_{s1}+P_{sv}{k}_{sv1} \\ {K}_{P_2}=P_{P_3}·k_{P_3-slave}=k_{v2}+P_s{k}_{s2}+P_{sv}{k}_{sv2} \end{gathered}$$

(7)

where the greater complexity of scattering (higher than second order) is not considered, and the propagation matrix related to the volume scattering mechanism is approximately equivalent to an identity matrix. The subscripts 1 and 2 represent the corresponding master image and slave image, and the subscripts $v$, $s$, and $sv$ represent the volume, surface, and dihedral scattering mechanisms, respectively. $P_s$ and $P_{sv}$ are the corresponding 3 × 3 propagation matrices of the surface and dihedral scattering. In the same way, $k_v$, $k_s$, and $k_{sv}$ are the scattering vectors of the volume, surface, and dihedral scattering mechanisms, respectively [32]. Moreover, the coherency matrix $T_6$ is formed by the 3-D coherent scattering vectors $K_{P_1}$ and $K_{P_2}$ for each end of the baseline, and is expressed as

$$\begin{gathered} T_6=\langle \left[\begin{array}{c}{K}_{P_1}\\ K_{P_2}\end{array}\right]·\left[\begin{array}{cc}{K}_{P_1}^H& K_{P_2}^H\end{array}\right]\rangle =\left[\begin{array}{cc}{T}_{11}& {\mathsf{\Omega }}_{12}\\ {\mathsf{\Omega }}_{21}& T_{22}\end{array}\right] \\ {T}_{11}=\langle K_{P_1}{K}_{P_1}^H\rangle ; T_{22}=\langle K_{P_2}{K}_{P_2}^H\rangle ;{\mathsf{\Omega }}_{12}=\langle K_{P_1}{K}_{P_2}^H\rangle \end{gathered}$$

(8)

With the combination of (7) and (8), the coherency matrices $T_{11}$ and ${\mathsf{\Omega }}_{12}$ can be decomposed into three components on the basis of the different scattering mechanisms, respectively, and are shown as

$$\begin{gathered} T_{11}=\langle K_{v_1}{K}_{v_1}^H\rangle +P_s\langle K_{s_1}{K}_{s_1}^H\rangle P_s^H+P_{sv}\langle K_{s{v}_1}{K}_{s{v}_1}^H\rangle P_{sv}^H \\ {\Omega }_{12}=\langle K_{v_1}{K}_{v_2}^H\rangle +P_s\langle K_{s_1}{K}_{s_2}^H\rangle P_s^H+P_{sv}\langle K_{s{v}_1}{K}_{s{v}_2}^H\rangle P_{sv}^H \end{gathered}$$

(9)

$$\begin{gathered} T_{11}=T_v+P_s{T}_s{P}_s^H+P_{sv}{T}_{sv}{P}_{sv}^H \\ {\Omega }_{12}={\Omega }_v+P_s{\Omega }_s{P}_s^H+P_{sv}{\Omega }_{sv}{P}_{sv}^H \end{gathered}$$

(10)

As a result, we can express the interferometric coherence of the single-polarization as

$$\tilde{\gamma }=\frac{{\omega }^H{\Omega }_{12}\omega }{\sqrt{{\omega }^H{T}_{11}\omega }·\sqrt{{\omega }^H{T}_{22}\omega }}\approx \frac{{\omega }^H{\Omega }_{12}\omega }{{\omega }^H{T}_{11}\omega }=\frac{{\omega }^H\left({\Omega }_v+P_s{\Omega }_s{P}_s^H+P_{sv}{\Omega }_{sv}{P}_{sv}^H\right)\omega }{{\omega }^H\left(T_v+P_s{T}_s{P}_s^H+P_{sv}{T}_{sv}{P}_{sv}^H\right)\omega }$$

(11)

where $\omega$ is the projection vector related to the polarization, $\tilde{\gamma }$ can be expressed as a function of the coherency matrices $T_i$ and ${\Omega }_i$ and the propagation matrix $P_i$ based on the approximate condition $\sqrt{{\omega }^H{T}_{11}\omega }·\sqrt{{\omega }^H{T}_{22}\omega }\approx {\omega }^H{T}_{11}\omega$. The subscript $i$ represents a kind of scattering mechanism. In this case, we can rewrite (11) based on the equivalent transformation, which is shown as [32]

$$\begin{gathered} {\omega }^H\left(P_i{\Omega }_i{P}_i^H\right)\omega =p_i{m}_{0i}{\tilde{\gamma }}_i \\ {\omega }^H\left(P_i{T}_i{P}_i^H\right)\omega =p_i{m}_{0i} \end{gathered}$$

(12)

where $m_0$ is the normalized radar cross section and $\tilde{\gamma }$ represents the coherence contribution. The propagation factor $p$ related to the propagation path can be described as a function of the azimuth observation angle ${\alpha }_{Az}$ and used to attenuate the corresponding contribution of each scattering mechanism $i$. Note that the sublook coherences are obtained from a pair of given full-resolution SAR images with a fixed polarimetric state. Therefore, the projection vector $\omega$ is constant for the different sublook coherences. As a result, the sublook interferometric coherences with different azimuth observation angles ${\alpha }_{Az}$ can be expressed as

$$\tilde{\gamma }\left({\alpha }_{Az}\right)=\frac{m_{0v}{\tilde{\gamma }}_v+p_s\left({\alpha }_{Az}\right)m_{0s}{\tilde{\gamma }}_s+p_{sv}\left({\alpha }_{Az}\right)m_{0sv}{\tilde{\gamma }}_{sv}}{m_{0v}+p_s\left({\alpha }_{Az}\right)m_{0s}+p_{sv}\left({\alpha }_{Az}\right)m_{0sv}}$$

(13)

In general, it is assumed that the surface and dihedral scattering mechanisms include equal coherence contributions, which is shown as

$${\tilde{\gamma }}_s={\tilde{\gamma }}_{sv}=e^{i{\phi }_0}$$

(14)

Therefore, the interferometric coherence related to a given polarimetric state can be rewritten, and is shown as

$$\tilde{\gamma }\left({\alpha }_{Az}\right)=e^{i{\phi }_0}\frac{m_{0v}{\tilde{\gamma }}_v{e}^{-i{\phi }_0}+p_s\left({\alpha }_{Az}\right)m_{0s}+p_{sv}\left({\alpha }_{Az}\right)m_{0sv}}{m_{0v}+p_s\left({\alpha }_{Az}\right)m_{0s}+p_{sv}\left({\alpha }_{Az}\right)m_{0sv}}=e^{i{\phi }_0}\frac{{\tilde{\gamma }}_{v0}+\mu \left({\alpha }_{Az}\right)}{1+\mu \left({\alpha }_{Az}\right)}$$

(15)

where ${\phi }_0$is the ground phase associated with the underlying topography, $\mu \left({\alpha }_{Az}\right)$ denotes the ground-to-volume scattering ratio, which is formulated as a function of the ground- (surface and dihedral) and volume- backscattering contribution, and is shown as

$$\mu \left({\alpha }_{Az}\right)=\frac{p_s\left({\alpha }_{Az}\right)m_{0s}+p_{sv}\left({\alpha }_{Az}\right)m_{0sv}}{m_{0v}}$$

(16)

In addition, ${\tilde{\gamma }}_{v0}$ represents the pure volume scattering coherence, which is modeled as a function of the forest’s biophysical parameters and the vertical wavenumber with an exponential backscattering profile, as shown in (17).

$$\begin{gathered} {\tilde{\gamma }}_{v0}\left(h_v,\sigma \right)={\tilde{\gamma }}_v{e}^{-i{\phi }_0}=\frac{{\int }_0^{h_v}f\left(z\right)e^{i{k}_{z}z}dz}{{\int }_0^{h_v}f\left(z\right)dz}=\frac{2\sigma \left(e^{\left(2\sigma h_v/cos\theta +i{k}_z{h}_v\right)}-1\right)}{\left(2\sigma +i{k}_{z}cos\theta \right)\left(e^{2\sigma h_v/cos\theta }-1\right)} \\ f\left(z\right)=e^{\frac{2\sigma h_v}{cos\theta }} \end{gathered}$$

(17)

where $h_v$ is the CHM and $\theta$ is the look angle of the master track. $\sigma$ denotes the extinction coefficient, which is used to describe the power attenuation of the SAR signals in the forest volume. $k_z$ is the vertical wavenumber impacted by the baseline parameters and topographic slope.

Finally, for the single-polarization InSAR data, the sublook interferometric coherence, which is related to the azimuth observation angle ${\alpha }_{Az}$, can be derived from the Equations (15)–(17) and is shown as

$$\tilde{\gamma }\left({\alpha }_{Az}\right)=e^{i{\phi }_0}\frac{{\tilde{\gamma }}_{v0}\left(h_v,\sigma \left({\alpha }_{Az}\right)\right)+\mu \left({\alpha }_{Az}\right)}{1+\mu \left({\alpha }_{Az}\right)}$$

(18)

where $\sigma$ is expressed as a function of the azimuth observation angle ${\alpha }_{Az}$. The reason is that ${\alpha }_{Az}$ results in the path difference of the sublook signals in the forest volume, and σ can be influenced by the azimuth observation angle ${\alpha }_{Az}$ in a forest volume with heterogeneous characteristics. Moreover, for the sublook coherences caused by the different azimuth observation angles, both the polarimetric state $\mathsf{\omega }$ and the variation interval of the azimuth observation angle $\Delta {\alpha }_{Az}$ are fixed. In this case, the corresponding ground-to-volume scattering ratio $\mu$ and extinction coefficient $\sigma$ can be influenced by the azimuth observation angle ${\alpha }_{Az}$. As a result, Equation (18), which has a similar expression to that of the OVoG model [34], is not invertible because of the underdetermined problem of unknown parameters.

To address this problem, the change in the extinction coefficient caused by the different azimuth observation angles within the volume scattering coherence ${\tilde{\gamma }}_{v0}$ can be ignored, since the difference in azimuth observation angles leads to a significant offset of the sublook coherences along the coherence line but a relatively small change in the pure volume coherence ($e^{i{\phi }_0}{\tilde{\gamma }}_{v0}$) in the complex plane, as has been verified in several studies [6,24]. Therefore, we can rewrite the sublook coherence as

$$\tilde{\gamma }\left({\alpha }_{Az}\right)=e^{i{\phi }_0}\frac{{\tilde{\gamma }}_{v0}\left(h_v,\sigma \right)+\mu \left({\alpha }_{Az}\right)}{1+\mu \left({\alpha }_{Az}\right)}=e^{i{\phi }_0}\left({\tilde{\gamma }}_{v0}+\frac{\mu \left({\alpha }_{Az}\right)}{1+\mu \left({\alpha }_{Az}\right)}\left(1-{\tilde{\gamma }}_{v0}\right)\right)$$

(19)

It is clear that Equation (19) shows the relationship between the sublook coherences and the forest’s biophysical parameters. According to the previous discussion, a full-resolution SAR image can be split into multiple sublook images by performing sublook decomposition, and then the corresponding sublook coherences can be acquired using interferometric processing (see Equation (2)). In this case, the number of observations made using the classical InSAR system is increased by changing the azimuth observation angle ${\alpha }_{Az}$, which makes it possible to quantitatively evaluate the unknown parameters based on the Equation (19). In the complex plane, the coherence points corresponding to different sublooks can be enveloped by a polygon called the coherence region [24]. The coherence region can be fitted as a line using the Line-fitting (LF) method, as in the case of the RVoG model. In such a case, it is possible to estimate the ground phase and the forest biophysical parameters in the complex circle.

Equation (19) describes the vegetation scattering model as a two-layer structure made up of a volume layer and a ground layer. It can be called the TF-RVoG model because of the use of the TF analysis. Note that the impacts of the different azimuth observation angles on the ground-to-volume scattering ratios are described by a pair of propagation factors $p_s\left({\alpha }_{Az}\right)$ and $p_{sv}\left({\alpha }_{Az}\right)$, as shown in Equation (16). The reason for this is that, with the difference in propagation paths caused by the varying azimuth observation angles, the sublook SAR signals undergo different attenuation in the heterogeneous forest volume, which results in a pair of varying propagation factors $p_s$ and $p_{sv}$. However, for the RVoG-based polarimetric observations, the ground-to-volume scattering ratio is determined by the different $m_{0s}\left(\omega \right)$, $m_{0sv}\left(\omega \right)$, and $m_{0v}\left(\omega \right)$. This indicates a significant difference in the impacts of the sublook and polarimetric coherences.

Moreover, an ideal case is provided to study in further detail the differences in the sublook and polarimetric InSAR observations, where the forest is a homogeneous volume with the same CHM and the underlying topography is a horizontal plane, as shown in Figure 2a. In this case, the pure volume coherence is a constant because of the fixed extinction coefficient of the homogeneous forest, and we can compare the coherence loci of the sublook and polarimetric observations in the complex unit circle. For the RVoG-based polarimetric coherences, the ground-to-volume scattering ratio is changed by the normalized radar cross section for the different polarimetric states. As a result, the coherence loci of the polarimetric observations strictly follow a straight line, as shown in Figure 2b. For the sublook coherences based on the TF-RVoG model, however, the propagation factors $p_s$ and $p_{sv}$ are invariable since the different sublook SAR signals undergo the same attenuation process in the ideal forest scene. Therefore, the corresponding coherence loci of the sublook observations present as isolated points with a constant ground-to-volume scattering ratio (Figure 2c).

2.2. RME-Induced Phase Error Correction

Due to the changes in the backscattering properties caused by the different propagation paths of the sublook SAR signals, the sublook coherences contain different interferometric phases. To accurately estimate the CHM from the sublook coherences, we must ensure that the change in the backscattering properties is the only factor leading to the interferometric phase difference of the sublook coherences. Therefore, the phase components that are induced by other factors varying with the different azimuth observation angles must be removed. For the airborne InSAR data used in this paper, the interferometric phase error caused by the residual motion error (RME) is the major influence factor [26].

With the impact of airflow turbulence, the airborne SAR sensor usually deviates from its predetermined track along the azimuth, which leads to a time-varying motion error. After interference processing, the motion error will be transformed into the error of the interferometric phase. In general, the motion error impacts of the different polarization coherences can be ignored when estimating the CHM by using the PolInSAR data. However, for the sublook InSAR coherences, a varying phase component will be induced by the different motion errors along the azimuth, which will have a negative impact on the CHM inversion. In other words, these motion errors are unabsorbable and must be removed for robust parameter estimation in the case of sublook InSAR. According to the existing studies, the motion parameters of aircraft can be record by the navigation system of the airborne SAR sensor, and most of the motion error can be removed through motion compensation processing [35,36,37]. However, there a fraction of RME will remain because of the limited precision of the current airborne navigation system.

To remove the RME-induced phase error of each sublook coherence, a parameterized model for the single-baseline case that can accurately describe the relationship between the errors of the time-varying baseline parameters and the RME-induced interferometric phase is used, and it can be expressed by Equation (20) [26]

$${\phi }_{diff}={\phi }_{RME}+{\phi }_{res-DEM}+{\phi }_{noi}=\left(\frac{\partial {\phi }_{res-flat}}{\partial B}+\frac{\partial {\phi }_{res-topo}}{\partial B}\right)·\Delta B+\left(\frac{\partial {\phi }_{res-flat}}{\partial \alpha }+\frac{\partial {\phi }_{res-topo}}{\partial \alpha }\right)·\Delta \alpha +{\phi }_0+\Delta \phi$$

(20)

where ${\phi }_{diff}$ is the differential interferometric phase, including three phase components: the RME-induced phase ${\phi }_{RME}$, the topographic error phase ${\phi }_{res-DEM}$ caused by discrepancies between the external digital elevation model (DEM) and the real DEM, and the noise phase ${\phi }_{noi}$. Furthermore, $\Delta B$ and $\Delta \alpha$ are the time-varying baseline parameter errors, which are related to the baseline length $B$ and the baseline tilt angle $\alpha$. Note that ${\phi }_{RME}$ can be divided into two parts, i.e., the residual flat-earth phase ${\phi }_{res-flat}$ and the residual topographic phase ${\phi }_{res-topo}$. ${\phi }_0$ represents the unwrapped phase of the reference point, and $\Delta \phi$ denotes the sum of the topographic error phase and the noise phase. Based on this model, three unknown parameters $\Delta B$, $\Delta \alpha$, and ${\phi }_0$ can be jointly estimated using a weighted least-squares (WLS) regression approach. As a result, we can acquire the corresponding RME-induced phase error and then deduct it from each sublook interferometric phase with the same processing.

2.3. Three-Stage Inversion Method Based on the TF-RVoG Model

The proposed TF-RVoG model has a similar expression to that of the RVoG model. However, for the sublook coherences based on the TF-RVoG model, it is difficult to estimate the CHM using the classical three-stage inversion method [15] for the following two reasons. Firstly, for the RVoG-based polarimetric observations, we can acquire the ground phase from the two intersection points between the coherence line and the unit circle in the complex unit circle using the a priori order of the polarimetric coherences along the coherence line. However, there is no a priori order of the sublook coherences to use in selecting the ground phase. Secondly, a measured coherence with a zero ground scattering component, i.e., the pure volume coherence, is needed to invert the CHM from a 2-D lookup table (LUT) with the CHM and the extinction coefficient. In the polarimetric cases, the coherence of HV or PDH polarization, which includes a lower ground scattering contribution in comparison with the other polarizations, is usually selected as the pure volume coherence. However, the pure volume coherence for the sublook observations is unknown.

To address these problems, a modified three-stage method that is suitable for the TF-RVoG model is proposed for the CHM inversion, and the detailed steps are as follows.

(1) Coherence Line Fitting: Based on Equation (19), the sublook coherences are impacted by the varying ground-to-volume scattering ratios caused by the different propagation factors. Hence, the coherence loci of the sublook InSAR observations form a finite line between the coherence of the ground phase and the pure volume coherence in the complex plane. Therefore, sublook coherences with different azimuth incidence angles can be fitted to a coherence line using the total least-squares (TLS) method.

(2) Ground Phase Selection: Two intersection points between the fitted coherence line and the complex unit circle are considered as candidate ground phase points, and we need to select the correct ground phase point from these two candidates without a priori knowledge. Due to the monotonic lowering of the sublook phases with the increase in the ground-to-volume scattering ratios, a robust criterion was derived in [38] that can be used to obtain the correct ground phase point in the sublook case. The corresponding decision rule is shown as

$$\begin{array}{ll}\mathrm{For}{k}_z<0:& \\ & \mathrm{if}\arg \left(\exp \left(i{\phi }_1\right)\exp {\left(i{\phi }_2\right)}^{*}\right)>0\mathrm{then}{\phi }_0={\phi }_1\\ & \mathrm{if}\arg \left(\exp \left(i{\phi }_1\right)\exp {\left(i{\phi }_2\right)}^{*}\right)<0\mathrm{then}{\phi }_0={\phi }_2\\ \mathrm{For}{k}_z>0:& \\ & \mathrm{if}\arg \left(\exp \left(i{\phi }_1\right)\exp {\left(i{\phi }_2\right)}^{*}\right)<0\mathrm{then}{\phi }_0={\phi }_1\\ & \mathrm{if}\arg \left(\exp \left(i{\phi }_1\right)\exp {\left(i{\phi }_2\right)}^{*}\right)>0\mathrm{then}{\phi }_0={\phi }_2\end{array}$$

(21)

where ${\phi }_0$ is the ground phase, and ${\phi }_1$ and ${\phi }_2$ are the two candidate ground phases. This criterion with a limited condition of $h_v\le \pi /k_z$ is effective for selecting the true ground phase and is free of the a priori order of the sublook coherences.

(3) Forest CHM Retrieval: To avoid the underdetermined problem in the CHM inversion, a coherence without the ground scattering contribution (i.e., $\mu \left({\alpha }_{Az}\right)=0$) needs to be selected from all the observed sublook coherences. Based on the opposite relationship between the interferometric phases of the sublook coherences and the corresponding ground-to-volume scattering ratios, the phase center height of each sublook coherence can be expressed as

$$d_i=\Vert \frac{arg\left(\gamma \left({\alpha }_{A{z}_i}\right)-{\phi }_0\right)}{k_z}\Vert$$

(22)

When the condition of $d_i=d_{max}$ is satisfied, the sublook coherence with the lowest ground scattering contribution can be acquired from the candidate sublook coherences. In this case, we can rewrite the (19) as

$$\tilde{\gamma }{\left({\alpha }_{Az}\right)}_{d_i=d_{max}}=e^{i{\phi }_0}{\tilde{\gamma }}_{v0}\left(h_v,\sigma \right)$$

(23)

As a result, a pair of the corresponding solutions $(h_v,\sigma$) can be retrieved in the form of the 2-D LUT.

2.4. 2-D Ambiguous Error Correction of the Pure Volume Coherence

For the sublook InSAR data, the CHM can be retrieved by the modified three-stage method without a priori knowledge. However, the accuracy of the CHM inversion is impacted by two major problems, i.e., the residual ground scattering error of the pure volume coherence and the temporal decorrelation error.

(1) Residual Ground Scattering Effect: Equation (22) shows that the sublook coherence with the lowest ground scattering contribution can be selected from all the sublook coherences, and it is regarded as the pure volume coherence. In the complex unit circle, the obtained pure volume coherence ${\gamma }_{v0}$ is located at the end of the coherence region farther from the ground phase point $e^{{\phi }_0}$, as shown in Figure 3a. However, due to the influence of the residual ground scattering error, the true pure volume coherence ${\gamma }_{v1}$ without considering the temporal decorrelation error should be located on a finite line between the ${\gamma }_{v0}$ and $e^{{\phi }_1}$, as shown in Figure 3b, which causes a one-dimensional (1-D) ambiguous error of the pure volume coherence along the coherence line.

(2) Temporal Decorrelation Effect: In an airborne repeat-pass interferogram, the temporal decorrelation error caused by repeated SAR acquisitions can result in an overestimated CHM [33,38,39]. During a short time interval (about 50 min in this paper), the wind-induced motion of the scatterers in the forest volume is the major factor causing the temporal decorrelation error of the airborne InSAR system. According to the existing studies, we can assume that the temporal decorrelation error only changes the amplitudes of the sublook coherences and has no effect on their phases [40]. In this case, the true pure volume coherence ${\gamma }_{v2}$ will be located on the extended segment of the finite line between ${\gamma }_{v1}$ and the center of the circle, forming the second 1-D ambiguous error (Figure 3c).

Therefore, in the complex unit circle, a 2-D ambiguous error in the pure volume coherence is caused by the joint influence of the residual ground scattering contribution and the temporal decorrelation, which needs to be removed to obtain a high-precision CHM.

To address these problems, firstly, a geometrical index $D.I$, which can interpret the relative location of the pure volume coherence on the coherence line, is expressed as in [25]

$$D.I=A.L/V.L$$

(24)

where $A.L$ denotes the ambiguous line length between $e^{{\phi }_1}$ and ${\gamma }_{v0}$, and $V.L$ represents the visible line length between ${\gamma }_{v0}$ and $e^{{\phi }_0}$, as shown in Figure 3a. Because of the monotonic lowering of the mean extinction with the increase in penetration depth, the pure volume coherence can be defined as a function of the extinction coefficient. In response, a relationship between the geometrical index $D.I$ and the extinction coefficient was proposed by performing a simulation experiment with an ideal pure volume coherence, as listed in

Table 1. The Limited Extinction Coefficient BASED ON the Geometrical Index.

Index Range	Extinction Coefficient Range (dB/m)
0–0.4	0.6–1
0.4–0.65	0.3–0.6
>0.65	0–0.3

[25], which has been proven to be effective for alleviating the impact of the residual ground scattering error.

Moreover, many methods have been proposed for removing the temporal decorrelation error [38,39,40]. A representative method that is based on the random volume over ground with volume temporal decorrelation (RVoG + VTD) model is to model the temporal decorrelation as a real-valued factor [40]. For the sublook case, the corresponding model can be expressed as

$$\tilde{\gamma }\left({\alpha }_{Az}\right)=e^{i{\phi }_0}\frac{{\alpha }_{vt}{\tilde{\gamma }}_{v0}\left(h_v,\sigma \right)+\mu \left({\alpha }_{Az}\right)}{1+\mu \left({\alpha }_{Az}\right)}$$

(25)

where ${\alpha }_{vt}$ represents the temporal decorrelation factor applied to the volume scattering coherence. In this case, the additional ${\alpha }_{vt}$, as an unknown parameter, leads to the underdetermined problem of the CHM inversion in the single-baseline InSAR case. This can be solved by fixing the extinction coefficient $\sigma$ or the temporal decorrelation factor ${\alpha }_{vt}$ [41]. However, both parameters are sensitive to the heterogeneity of the forest’s vertical structure (e.g., forest CHM, forest density, and forest species). As a result, a non-negligible error could be caused by the fixed extinction coefficient or temporal decorrelation for the CHM inversion. In this paper, the relationship between the CHM and temporal decorrelation is studied in light of the spatial heterogeneity of forest scenes. Due to the short space of time in which the repeated SAR acquisitions were made in the airborne case, most of the temporal decorrelation is attributed to the wind-induced motion of scatterers. Note, however, that the influence of wind can lead to a greater temporal decorrelation error for the high forest, compared with the low forest. The reason for this is that the high forest with its larger leaves and longer twigs is more susceptible to the wind in comparison with the low forest. Therefore, a linear relationship can be fitted to describe the interaction of temporal decorrelation factor with the CHM (i.e., ${\alpha }_{vt}=a{h}_v+b$). As a result, Equation (25) can be rewritten and called the random volume over ground with mutative temporal decorrelation (RVoG + MTD) model, and shown as

$$\tilde{\gamma }\left({\alpha }_{Az}\right)=e^{i{\phi }_0}\frac{\left(a{h}_v+b\right){\tilde{\gamma }}_{v0}\left(h_v,\sigma \right)+\mu \left({\alpha }_{Az}\right)}{1+\mu \left({\alpha }_{Az}\right)}$$

(26)

For the RVoG + VTD model, we calculated the unknown extinction coefficient (or temporal decorrelation) by using a small part of the a priori CHM, and then fixed it with an average. However, for the RVoG + MTD model in Equation (26), there are three unknown parameters ($a$, $b$, and $\sigma$) that need to be retrieved. In order to avoid the underdetermined problem, the inversion process is divided into three stages. Firstly, the temporal decorrelation factor is derived using the a priori CHM. Moreover, the unknown parameters $a$ and $b$ can be estimated by fitting the linear relationship between the temporal decorrelation factor and the a priori CHM. Finally, we inverted the corresponding solution pairs ($h_v,\sigma$) using the 2-D LUT.

In conclusion, with the limitation of the empirical relationship between the $D.I$ and the extinction coefficient (see Table 1), the proposed RVoG + MTD model is used to alleviate the impact of the 2-D ambiguous error of pure volume coherence and to improve the accuracy of the CHM retrieval.

3. Experiments and Results

In this study, the performance of the proposed approach to CHM inversion (Figure 4) was tested with the L-band full-polarization SAR data acquired by the airborne E-SAR system of the German Aerospace Center (DLR) during the BioSAR 2008 campaign.

3.1. Experimental Data

The test site is a forested area in the Vindeln municipality, located in the Krycklan catchment of Northern Sweden ($64°{10}^{\prime }$N and $20°{01}^{\prime }$E). In this study area, the main forest types are coniferous tree species (Norway spruce and Scots pine), which cover about 70% of the total area. The rest is broad-leaved forest. The forest height varies from 0 to 35 m with a mean forest height of 18 m and a mean biomass level of 90 t/ha, and the underlying ground height ranges from 150 to 380 m above mean sea level with a moderate slope ranging from $-15°$ to $15°$, as shown in Figure 5a. For the L-band fully polarimetric SAR data, one interferometric pair with a spatial baseline of 12 m was used to retrieve the CHM, and the time interval of the repeated SAR acquisitions was about 50 min. The SAR images have the pixel spacing of 0.74 m in azimuth and 1.5 m in slant range; Figure 5b shows the corresponding Pauli-basis image. In addition, the airborne LiDAR CHM, collected by the Swedish Defense Research Agency (FOI) in late summer 2008, is used to evaluate the accuracy of the InSAR-derived CHM and is shown in Figure 5c.

3.2. Processing of the Sublook InSAR Data

In order to perform the single-polarization sublook InSAR experiment, the HH polarization was selected from the observed full-polarization SAR data for the following reasons. Firstly, the co-polarization HH or VV is usually used in the single-polarization InSAR system, due to the higher coherence in comparison with the cross-polarization (HV and VH). In addition, the coherence of the HH polarization, which includes the lower volume scattering contribution, is higher than the coherence of the VV polarization.

According to the Doppler bandwidth criterion, the sublook images with a wide bandwidth have a high resolution but with a small coherence dispersion. In contrast, the sublook images with a narrow bandwidth result in a low resolution, but the coherence dispersion is larger. In this case, a suitable bandwidth needs to be selected for a robust parameter inversion. According to the analysis of the Doppler Bandwidth described in [24], in this paper, the master and slave images of the HH polarization were decomposed into five sublook images with 50% of the total Doppler bandwidth, respectively. Meanwhile, for a fair comparison between the polarization images and the sublook images, an equivalent number of looks is necessary. For this purpose, the polarization coherences were calculated over an 11 × 11 sliding window, while the sublook coherences were estimated with a 15 × 15 sliding window [30].

After interferometric processing, the phase error of each sublook coherence caused by the RME can be removed using the time-varying baseline parameter errors model, as shown in Equation (20). The interferometric phases of the five sublook coherences need to be unwrapped before removing the RME-induced errors. Figure 6 shows the results of the phase unwrapping of the first sublook coherence. The coherence and phase maps of the first sublook observation were shown in Figure 6a and Figure 6b, respectively. The minimum cost flow (MCF) method was used to unwrap the interferometric phase [42], and the unwrapped phase of the first sublook coherence is shown in Figure 6c. For the study area shown in the red polygons, a low coherence ($\gamma <0.45$) resulted in unsmooth phases during the unwrapping process. This problem can be alleviated by using the difference phases from the different pairs of sublook coherences, and the optimized result of Figure 6c is shown in Figure 6d.

In addition, the parameterized error model for the time-varying baseline parameters can be used to estimate the RME-induced phase error of each sublook coherence. As a result, the phase errors of the first and the fifth sublook coherences caused by the RME are shown in Figure 7a and Figure 7b, respectively. To present the impacts of the RME-induced phase errors clearly, we have studied the difference maps between the first and fifth sublook phases (i.e., ${\phi }_5-{\phi }_1$) before and after the RME-induced phase errors were removed, as shown in Figure 7c,d. For the difference map shown in Figure 7c, it is clear that there are some phase fluctuations along the azimuth. The major reason for this is that the phase errors of the sublook coherences induced by the RME are different. Therefore, the corresponding RME-induced phase error of each sublook coherence should be removed, and the improved result of the difference phase is shown in Figure 7d.

3.3. Results and Analyses

The obtained CHM maps based on the different methods from the PolInSAR interferometric pair (with a spatial baseline of 12 m and a temporal baseline of about 50 min) are shown in Figure 8a–f scaled from 0 to 35 m. Moreover, in order to obtain reliable forest CHM results from the single baseline acquisitions, a general compromise is to mask the areas which have too large or small a value for the vertical wavenumber [38], i.e., the areas with $k_z<0.06$ rad/m and $k_z>0.2$ rad/m in the present study. In the available area, 181 out of the 256 cross-validation stands are superimposed on the CHM map acquired by PolInSAR (Figure 8a). To more clearly present the contribution of the proposed models and methods, the CHM inversion is described as three stages, and more details are given as follows.

(1) CHM Inversion based on the TF-RVoG model: For a fair comparison of the different methods, the classical three-stage inversion method was used to retrieve the CHM from the PolInSAR data [15], which can be regarded as a reference to verify the performance of the proposed approach. In this case, however, we must be aware that the polarization selection has a great impact on coherence line fitting and parameter retrieval for the RVoG model. For the linear polarization coherences, their high similarity in the complex plane may lead to fitting errors of the coherence line. In contrast, the optimized coherences can improve the linearity of the coherence sets because of the greater differences in the coherences along the coherence line. Therefore, five optimized coherences, namely the phase diversity (PD) coherences PDH and PDL [43], and the optimum coherences Opt1, Opt2, and Opt3 [44], are used to fit the coherence line. As a result, the RVoG-derived CHM is shown as Figure 8a. In addition, according to the modified three-stage inversion method in this study, the corresponding CHM can be derived using the proposed TF-RVoG model from five sublook InSAR coherences with the same bandwidths (50% of the total Doppler bandwidth) in HH polarization and is shown as Figure 8b.

Compared with the CHM mapped by the LiDAR, as shown in Figure 8g, it is clear that the CHMs of both Figure 8a and Figure 8b show obvious overestimation. The major reason for this is that the joint influence of the residual ground scattering contribution of the pure volume coherence (i.e., the coherence of PDH polarization in the PolInSAR case or the coherence of $d_i=d_{max}$ in the sublook InSAR case) and the temporal decorrelation. The 181 cross-validation plots of the estimated CHM against the LiDAR-derived CHM are displayed in Figure 9, and the corresponding root-mean-square error (RMSE) and determination coefficient ($R^2$) can be calculated. Firstly, the CHM inversion results derived by the PolInSAR RVoG method and the InSAR TF-RVoG method reach a determination coefficient of 0.62/0.64 with an RMSE of 6.55/5.61 m, as shown in Figure 9a and Figure 9b, respectively. This indicates that we can retrieve the CHM from the single-baseline InSAR data without a priori knowledge, and the CHM result obtained using the InSAR TF-RVoG method shows an obvious improvement of 14.3% in comparison with the RVoG-derived result from the PolInSAR data. Secondly, the high forest (15 to 35 m) has a greater overestimation than the low forest (0 to 15 m). The reason for this is likely that the higher forest has a weaker resistance to the errors in the pure volume coherence. Finally, the consistency of the inversion results obtained by the above-mentioned two methods is underlined by a direct comparison of the CHM estimates with an $R^2$ of 0.87 and an RMSE of 2.52 m, as shown in Figure 9c.

In order to further verify the effectiveness of the proposed method, we studies the coherence differences between the PolInSAR data and the sublook InSAR data in two major aspects, i.e., the coherence line fitting of the measured coherences and the residual ground scattering contribution of the selected pure volume coherence. As a result, the relevant indexes, including CLI, prod, and PCHV (i.e., phase center height of the volume-only coherence), are used for the coherence analysis in the complex plane. Firstly, the CLI can be used to describe the linear consistency of the five measured coherences, and is expressed as Equation (27) [6]

$$CLI=\frac{{\sum }_{k=1}^n{d}_k}{L}$$

(27)

where $d$ denotes the vertical distance from a certain coherence point to the fitted coherence line, and $L$ represents the distance between the two vertical projection points on the coherence line, which is induced by the two ends of the coherent region. $n$ is the number of coherence points, and $n=5$ in this study. Moreover, the prod is an index used to describe the maximum coherence difference and the mean magnitude of the coherence region, and it can be formulated as Equation (28) [17]

$$prod=\left|{\gamma }_{high}-{\gamma }_{low}\right|·\left|{\gamma }_{high}+{\gamma }_{low}\right|$$

(28)

where $\left|{\gamma }_{high}-{\gamma }_{low}\right|$ and $\left|{\gamma }_{high}+{\gamma }_{low}\right|$ represent the length of the major axis line of the coherence region and two times the magnitude of the center of the coherence region, respectively. In general, a fitted coherence line with both a low value for the CLI and a high value for the prod helps acquire a high-accuracy ground phase. Finally, the PCHV is used to express the influence of the residual ground scattering error of the pure volume coherence, and can be expressed as Equation (22). For the pure volume coherence, a higher PCHV corresponds to a lower error induced by the residual ground scattering contribution.

In this study, these indexes, such as CLI, prod, and PCHV, can be estimated from the measured coherences of the PolInSAR and the sublook InSAR, respectively. Figure 10a shows the CLI histograms obtained by the two methods mentioned. Clearly, the CLI result derived using the InSAR TF-RVoG method has a smaller mean than the PolInSAR RVoG-derived result. In other words, the measured sublook coherences present a higher linear correlation than the optimized polarization coherences. Figure 10b shows a big difference in the distribution frequency of the prod results between the two methods, but the close means of 0.358/0.303 are calculated separately from the coherences of the PolInSAR and the sublook InSAR. Furthermore, the comparison of the PCHV results, as shown in Figure 10c, indicates that it is reasonable that the coherence with $d_i=d_{max}$ (i.e., SVC) is used as the pure volume coherence in the sublook InSAR case, where it plays a similar role to that of the PDH coherence in the PolInSAR case. In conclusion, the proposed sublook TF-RVoG method can be used to retrieve the CHM with an obvious improvement in comparison with the existing PolInSAR RVoG method, and it is free from the multi-polarization and the observation efficiency reduction.

(2) Correction of the residual ground scattering error: The 1-D ambiguous error of the pure volume coherence caused by the residual ground scattering contribution (Figure 3b) can be alleviated using a limited extinction coefficient derived from the empirical relationship between the geometrical index $D.I$ and the extinction coefficient (Table 1). Figure 8b,c show the CHMs derived using the proposed InSAR TF-RVoG approach before and after the empirical relationship is used, respectively. On this basis, Figure 11 shows the cross-validation results for the InSAR-derived CHM vs. the LiDAR-derived CHM. It is clear that, with the limitation of the empirical relationship, the determination coefficient improves slightly with an $R^2$ of 0.71, whereas the RMSE reduces to 4.01 m.

To investigate the detailed differences in the CHM inversion results before and after the empirical relationship is used, the coherence effect caused by the limited extinction coefficient is expressed in the complex unit circle (Figure 12). Firstly, Figure 12a shows how to search for the pure volume coherence based on the 1-D ambiguous loci between ${\gamma }_v$ and $e^{i{\phi }_1}$ along the coherence line, where ${\gamma }_v$ denotes the pure volume coherence regardless of the residual ground scattering error. The red curve represents the corresponding loci of the coherence ${\gamma }_i$ with a fixed extinction coefficient σ and the varying CHM $h_v$. We can obtain a unique solution for $h_v$ for the condition ${\gamma }_i={\gamma }_v$. In this case, however, the $h_v$ is usually overestimated because of the residual ground scattering effect, and the true pure volume coherence should be located to the left of ${\gamma }_v$ in the coherence line. In response, the $D.I$ is calculated, and the $\sigma$ can be updated using the empirical relationship of $D.I$ and $\sigma$. As a result, the blue curve shows the coherence loci with a fixed $\sigma$ after updating, and ${\gamma }_v^{\text{'}}$ represents the improved pure volume coherence, which is the intersection point between the blue curve and the coherence line. Moreover, Figure 12b shows the coherence regions divided by the empirical relationship of $D.I$ and $\sigma$, where the green dots represent the coherence loci obtained from the varying $D.I$ from 0 to 0.4. In this case, the corresponding σ ranges from 1 dB/m to 0.6 dB/m. In a similar way, the blue dots represent the coherence loci obtained from the varying $D.I$ from 0.4 to 0.65, and the red dots the coherence loci obtained from the varying $D.I$ greater than 0.65. The corresponding $\sigma$ range from 0.6 dB/m to 0.3 dB/m, and from 0.3 dB/m to 0, respectively.

The CHM result derived using the improved method with the limited $\sigma$, as shown in Figure 11b, is still overestimated for the high CHM. The reasons are likely that, on the one hand, the temporal decorrelation is a non-negligible factor in CHM overestimation when the L-band InSAR data are used; on the other hand, the ground scattering error is not completely removed, due to the rough ranges of the $D.I$ (i.e., $0<D.I\le 0.4$, $0.4<D.I\le 0.65$, and $D.I>0.65$). In response, some representative points from the experimental data are presented for detailed description, as listed in

Table 2. Comparison of the CHMs before and after using the limited extinction coefficient.

$\mathit{D}.\mathit{I}$	InSAR CHM (m)	InSAR CHM with a Limited $\mathit{\sigma }$	$\mathbf{Distance}\mathbf{between}{\mathit{\gamma }}_{\mathit{v}}$$\mathbf{and}{\mathit{\gamma }}_{\mathit{v}}^{\prime }$
0.4	14.5	9.5	0.126
0.32	21	14.5	0.267
0.43	20.5	14	0.198
0.71	35	35	0
0.22	17	12.5	0.055
0.23	22.5	15	0.054
2.02	6.5	6.5	0
0.82	32	32	0

. It is clear that, for the case of $D.I\le 0.65$, the corresponding range of the limited $\sigma$ leads to the selection of the pure volume coherence ${\gamma }_v^{\text{'}}$ with the lower ground scattering contribution (Figure 12a) and the lower $h_v$. However, for the cases of $D.I>0.65$, the coherences can be shown as the red dots in the Figure 12b. In this case, ${\gamma }_v^{\text{'}}={\gamma }_v$, and the $h_v$ remains unchanged. In other words, the restricted relationship of $D.I$ and $\sigma$ is ineffective for removing the residual ground scattering error of the pure volume coherence in the case of $D.I>0.65$.

(3) Correction of the temporal decorrelation error: The wind-induced temporal decorrelation error can be removed by using the RVoG + MTD model related to a linear relationship of ${\alpha }_{vt}=a{h}_v+b$, as shown in Equation (26). Firstly, 12 validation stands obtained using LiDAR are used as the a priori CHM [see Figure 13a]; then, the corresponding temporal decorrelation factor ${\alpha }_{vt}$ can be calculated from the RVoG + VTD model [see Equation (25)] with a 2-D LUT of ${\alpha }_{vt}$and σ, as shown in Figure 13b. Comparing the results between Figure 13a and the Figure 13b, it is clear that increasing the CHM results in a lowering of the temporal decorrelation factor. This indicates that, for the high forest, greater amounts of temporal decorrelation error are induced. In contrast, the areas with a smaller error of the temporal decorrelation include the low forest. Secondly, for the 12 selected validation stands, the linear relationship of the CHM and the temporal decorrelation factor (i.e., ${\alpha }_{vt}=a{h}_v+b$) can be fitted in the limited ranges of the $D.I$ (i.e., $0<D.I\le 0.4$, $0.4<D.I\le 0.65$, and $D.I>0.65$), and the corresponding results are shown in Figure 14a–c, respectively. As a result, the temporal decorrelation factor, as a function of CHM, can be expressed as

$$\{\begin{array}{cc}0<D.I\le 0.4,& {\alpha }_{vt}=-0.022h_v+1.006;\\ 0.4<D.I\le 0.65,& {\alpha }_{vt}=-0.035h_v+1.257;\\ D.I>0.65,& {\alpha }_{vt}=-0.031h_v+1.205;\end{array}$$

(29)

Therefore, we can replace the temporal decorrelation factor ${\alpha }_{vt}$ with the equivalent expression $h_v$ (Equation (29)). In this case, for the study area, the unknown solution pairs ($h_v,\sigma$) can be retrieved by using the RVoG + MTD model.

To verify the performance of the proposed approach, two existing approaches based on the RVoG + VTD model, namely fixing $\sigma$ or fixing ${\alpha }_{vt}$ as a constant value, are used as references. It must be emphasized that a small amount of a priori CHM is necessary for both of the methods to obtain the fixed $\sigma$ or ${\alpha }_{vt}$. Therefore, the 12 validation stands, as shown in Figure 13a, are used as the a priori CHM for a fair comparison of the different methods. Moreover, for each method, the pure volume coherence error induced by the residual ground scattering effect has been alleviated by using the empirical relationship between $D.I$ and $\sigma$. Subsequently, we compared the impacts of the temporal decorrelation errors for different methods, including the RVoG + VTD model with a fixed $\sigma$, the RVoG + VTD model with a fixed ${\alpha }_{vt}$, and the RVoG + MTD model with a varying ${\alpha }_{vt}$.

Figure 15 shows the cross-validation results for the InSAR-derived CHM versus the LIDAR-derived CHM. Firstly, the resulting scatter plot derived using the RVoG + VTD model with a fixed σ is shown in Figure 15a. It is clear that the estimated CHM exhibits underestimation, and that the RMSE is 4.54 m and the $R^2$ is 0.51. The reason is that, with the limitation of the relationship between $D.I$ and $\sigma$, $\sigma$ should be divided into three different ranges (i.e., 0 < σ < 0.3, 0.3 ≤ σ < 0.6, and 0.6 ≤ σ < 1), and is fixed with an average, which leads to the overestimation of temporal decorrelation error and the underestimation of CHM. Secondly, a better CHM result is acquired by the RVoG + VTD with a fixed ${\alpha }_{vt}$, as shown in Figure 8e, and it is assessed in Figure 15b. The corresponding RMSE is 3.78 m with an $R^2$ of 0.71. In this case, for the high forest ranges from 20 m to 35 m, the InSAR-derived CHM presents an overestimation, but the low forest (0 to 20 m) shows obvious underestimation. This is because the taller trees yield greater amounts of temporal decorrelation, and the fixed ${\alpha }_{vt}$ is derived by an average of ${\alpha }_{vt}$, causing the underestimation and the overestimation of the temporal decorrelation error in the high and low forest, respectively. Finally, in the RVoG + MTD case, the spatial heterogeneity of the forest height is considered in the proposed approach by using a varying ${\alpha }_{vt}$, as shown in Equation (28), which helps to accurately remove the temporal decorrelation error. As a result, the CHM is estimated with an $R^2$ of 0.77 and an RMSE of 2.54 m [see Figure 15c], which represents an average improvement of 32.8% in comparison to the CHM obtained by the RVoG + VTD with a fixed ${\alpha }_{vt}$.

4. Discussion

The cross-validation results for the CHMs derived using the different methods, as shown in Figure 9, Figure 11 and Figure 15, have demonstrated the validity of the proposed model and method. In other words, the proposed CHM inversion approach was capable of estimating the CHM from the single-baseline InSAR data and improved both the accuracy and the robustness of the forest height retrieval. However, some potential problems with the proposed approach need to be further discussed and studied in order to improve the performance of the proposed approach.

4.1. The Limitation of the External DEM in the RME Correction

In order to obtain a robust CHM from the sublook InSAR data, the RME-induced phase error of each sublook coherence needs to be removed using the RME correction method. The weighted phase curvature autofocus method is free from problems with low coherence and can estimate the RME for an individual SAR image, but it requires high signal-to-noise ratio (SNR) targets in the SAR image [45]. The interferometric phase calibration approach can calibrate the phase errors by detecting stable point-like targets, but multi-baseline InSAR data are required [46]. In contrast, the time-varying baseline parameter errors model used in this paper is not only applicable to an interferogram characterized by low-coherence regions, but can also work on single-baseline InSAR data without requiring high-SNR targets. In addition, the topographic phase, which is important in detecting the RME-induced phase from the airborne repeat-pass interferometric phase, is related to the external DEM. Therefore, the performance of the time-varying baseline parameter errors model may be impacted by the error of the external DEM. Note that, in this study, the LiDAR-derived DEM with high resolution and high accuracy is used as the external DEM. However, for other study areas, it is difficult to find a DEM that matches the airborne InSAR data. The openly released spaceborne DEM products, such as SRTM DEM, ASTER DEM, and TanDEM-X DEM, can also play the role of the external DEM, but their effectiveness is limited because of their low resolution and precision [47]. Fortunately, the accuracy of the estimated CHM is not sensitive to the impact of the external DEM error in the proposed approach.

4.2. The Impact of the Rough Ranges of the $D.I$ in the Ground Scattering Error Removal

In this study, with the limitation of the relationship between the geometrical index $D.I$ and the extinction coefficient, the coherence regions can be divided into three parts, as shown in Figure 12b. In this case, an extinction coefficient with a higher accuracy can be acquired, which helps alleviate the residual ground scattering error of the pure volume coherence. However, for the representative points from the study area listed in Table 2, if $D.I>0.65$, the pure volume coherence remains unchanged before and after the empirical relationship is used, and the corresponding accuracy of the CHM has not been improved. This means that the empirical relationship of $D.I$ and the extinction coefficient only narrows the search range for the extinction coefficient, but does not provide its exact value. Therefore, the ground scattering errors are not completely removed, due to the rough ranges of the $D.I$. In order to improve the accuracy of the CHM inversion, a more precise relationship between the $D.I$ and the extinction coefficient needs to be further explored in future work.

5. Conclusions

In this study, the path difference of the sublook SAR signals, which is the major factor causing the backscattering variations, was considered in the scattering matrix. On this basis, a TF-RVoG model has been developed from the existing RVoG model and can be used to describe the relationship between the sublook coherence and the forest biophysical parameters. The corresponding three-stage method has been proposed to allow us to estimate CHM from the single-baseline sublook InSAR data. The inversion scheme was validated using L-band HH-polarization InSAR data. The RMSE of the validated plots based on the proposed method was 5.61 m, representing a slight improvement of 14.3% in comparison with the classical three-stage inversion result from the PolInSAR data. This indicates the CHM can be retrieved from the single-baseline InSAR data without a small amount of a priori CHM.

In addition, for a higher accuracy of the CHM, the impact of the residual ground scattering error of pure volume coherence was alleviated by using a limited extinction coefficient derived from the empirical relationship of the geometrical index $D.I$ and the extinction coefficient. As a result, the RMSE is reduced to 4.01 m, and the determination coefficient improves, with an $R^2$ of 0.71. To reduce the impact of the temporal decorrelation error in the CHM inversion, the relationship of the CHM and the temporal decorrelation is studied in view of the spatial heterogeneity of the CHM. Based on the negative correlation between the CHM and the temporal decorrelation, we describe the temporal decorrelation factor as a function of the CHM using linear fitting, and propose the RVoG + MTD model with a varying temporal decorrelation factor. Finally, a CHM is retrieved with an RMSE of 2.54 m, which represents an average improvement of 32.8% in comparison to the case with a fixed temporal decorrelation factor. The results indicate that, in combination with the relationship of the $D.I$ and the extinction coefficient, the proposed RVoG + MTD method is an effective way to alleviate the 2-D ambiguous error of the pure volume coherence induced by the joint impact of the residual ground scattering and the temporal decorrelation.

The performance of the proposed approach is hindered by the errors of the external DEM and the rough division of the $D.I$, which reduce the accuracy of the CHM. Moreover, as in existing methods, a small amount of a priori CHM is needed for the RVoG + MTD method to remove the error of the temporal decorrelation, which limits the scope of application of the proposed method. To overcome these problems, we will focus on how to use the proposed method with multi-baseline InSAR data, which has the potential to further improve the accuracy of the CHM owing to its greater number of available observations.