Forest Height Estimation Using a Single-Pass Airborne L-Band Polarimetric and Interferometric SAR System and Tomographic Techniques

Abstract

This paper addresses forest height estimation for boreal forests at the test site of Edson in Alberta, Canada, using dual-baseline PolInSAR dataset measured by Intermap’s single-pass system. This particular dataset is acquired by using both ping-pong and non-ping-pong modes, which permit forming a dual-baseline TomoSAR configuration, i.e., an extreme configuration for tomographic processing. A tomographic approach, based on polarimetric Capon and MUSIC estimators, is proposed to estimate the elevation of tree top and of underlying ground, and hence forest height is estimated. The resulting forest DTM and DSM over the test site are validated against LiDAR-derived estimates, demonstrating the undeniable capability of the single-pass L-band PolInSAR system for forest monitoring.

Short Biography of Author

Dr. Qiaoping Zhang is a geomatics scientist with more than 20 years of research and development experience in geographic information system, photogrammetry, and remote sensing. He received his first doctoral degree in Photogrammetry and Remote Sensing fromWuhan University (China) in 2002 and a second doctoral degree in Geomatics Engineering from University of Calgary (Canada) in 2006. Since 2006, he has been one of the leading researchers at Intermap. He has co-authored more than 40 journal research papers and conference publications. He was one of the ASPRS GeoEye Foundation Award recipients in 2006 and he held Alberta Innovates Fund–Research and Development Associate from 2006 to 2008. He is an active member of ASPRS and the Association of Professional Engineers, Geologists and Geophysicists of Alberta.

1. Introduction

Forests play a crucial role in Earth’s carbon cycle by absorbing carbon from the atmosphere and storing it in its biomass. Forest characterization and monitoring are hence very important for tracking climate change and quantifying the global carbon cycle in the form of above-ground-biomass mapping. SAR remote sensing techniques that guarantee under-foliage penetration are able to provide reliable estimation of forest parameters on a global scale [1], with a suitable temporal and spatial resolution. Among the SAR-extracted forest parameters, forest height is a key descriptor, commonly used for the estimation of forest volume and above-ground biomass. Polarimetric SAR Interferometry (PolInSAR), together with an appropriate inversion model such as the Random Volume Over Ground (RVOG) [2,3,4], has demonstrated a significant success for the recovery of canopy height and of the underlying ground elevation. Using Multi-baseline PolInSAR techniques, SAR tomography (TomoSAR) offers an efficient tool to carry out accurate 3D reconstruction of forested areas [5,6,7,8,9,10,11], including forest vertical structure and physical features. Forest height estimation has been successfully achieved via some well-known repeat-pass L-band PolInSAR systems in the literature. For example, the established single-baseline PolInSAR inversion techniques have been used to estimate tree heights over temperate forests, using L-band PolInSAR data acquired by DLR’s F-SAR system [12], as well as over tropical forests using L-band PolInSAR data acquired by E-SAR sensor during the Indrex II campaign [13]. TomoSAR techniques have also been applied to a six-pass dataset acquired by E-SAR sensor during BioSAR campaign over boreal forests [14].

In 2007 and 2008, Intermap developed an experimental single-pass L-Band PolInSAR system operated at a relatively lower altitude (∼1000 m) with a physical baseline around 3.5 m. The design philosophy was driven by the desire to demonstrate, relatively quickly and inexpensively, the capabilities and issues associated with tree height and Digital Terrain Model (DTM) estimation using PolInSAR techniques. The single-pass configuration permits assessing performance in forest monitoring in the absence of temporal decorrelation and residual motion issues. Mercer et al. [15] reported RVOG inversion results from single-baseline datasets, regarding DTM extraction and then forest height estimation over areas largely populated by pine stands, with heights ranging up to 30 m. Those results suggest that tree height accuracy (relative to LiDAR-derived heights) in sampled subsets was better than 10% of tree height, while the DTM showed biases that were typically 2–3 m above the true ground elevation. The RVOG-based PolInSAR retrieval techniques also met some limitations at the near range of the acquired images, which reduced the useful swath of this system for forest monitoring.

Regarding this PolInSAR system, a single transmitter and receiver chain provides fully polarimetric data through a switching network that allowed pulses to be recorded in both ping-pong and non-ping-pong modes simultaneously. Three independent antenna phase center positions can hence be obtained in a single flight, leading to a dual-baseline configuration, which provides the feasibility of applying SAR tomography to forest height estimation. This paper aims to evaluate the performance of forest height estimation using this dual-baseline TomoSAR configuration. Classical tomographic methods, such as Capon and MUSIC techniques, may meet some limitations for 3D reconstruction of forested areas, when dealing with this TomoSAR configuration whose vertical resolution is coarser than the observed tree heights. To overcome these limitations, the Digital Surface Model (DSM) of forests and the underlying ground elevation (DTM) are both estimated by the proposed tomographic techniques, and then validated against LiDAR data.

This paper is organized as follows. A detailed presentation of our test site and acquisition system is provided in Section 2. Tomographic processing for estimation of forest DTM and DSM is described in Section 3, including the limitations of two classical methods with respect to the present dual-baseline TomoSAR configuration as well as the proposed solution. Section 4 illustrates a solution of model order selection used for the proposed tomographic approach, using forest/non-forest mapping. Section 5 depicts the DTM, DSM and tree heights estimated by the proposed tomographic techniques, as well as those validated against the LiDAR measurements.

2. Presentation of Study Data and Acquisition System

2.1. Test Site Presentation

The test site is a forested area near the town of Edson in central Alberta, Canada (Figure 1a) [15]. This area is well suited for experiments related to forest height estimation, as it consists of a pattern of alternating forests and clear-cut patches. The forest over the test site is mainly composed of lodge pole pine trees with heights ranging from 10 to 30 m and with relatively homogeneous growth in the forested patches. Figure 1b shows ground photos over different patches of forests as well as short vegetation areas at the test site. The scene topography was measured using a LiDAR device (courtesy Terrapoint Canada) and the corresponding datasets were processed, in order to compute DSM and DTM estimates. In particular, the tree top height is derived from the envelope of the measured data in Figure 2, using the highest point within a 5.6 m searching radius, which results in an indicator similar to the H100 metric [13].

2.2. Acquisition System

The radar system considered in this study consists of an original quad-polarization L-band (wavelength $\lambda =0.2262$ m) radio-frequency hardware, implemented on Intermap’s TopoSAR platform. Supported by a rigid beam structure, two log-periodic antennas are mounted at the ends of a baseline $B=3.5$ m. Basic interferometric operation is carried out in ping-pong mode with both antennas transmitting and receiving separately, which in combination with polarimetric operation adds up to 8 individual channels (HH, HV, VH, and VV for each antenna). To add a second interferometric baseline, the system can be additionally operated in bistatic mode, i.e., one antenna transmitting and the other one receiving, which essentially halves the physical baseline if combined with one of the ping-pong channel sets. With this configuration, the interferometer dynamic range with respect to vegetation penetration is being increased significantly. Regarding the main system parameters, the most notable is the high geometric resolution of 1.1 m in range and 0.5 m in azimuth that allows efficient incoherent averaging and speckle filtering, while maintaining a decent product resolution. The comparatively low flight altitude (1000 m) and narrow swath width is a consequence of the low-cost experimental setup, using a short physical baseline, low power, and low-gain antennas.

2.3. Single-Pass Dual-Baseline Tomographic Configuration

TomoSAR, with an additional aperture in elevation formed by multiple acquisitions, is able to reconstruct the 3D structures of semi-transparent media, and in particular forests. Tomographic techniques can provide very good estimates of tree heights and underlying ground topography using repeat-pass Fully Polarimetric (FP) datasets [8,16].

The acquisition system used in this study offers the possibility to acquire Multi-Baseline Polarimetric Interferometric SAR (MB-PolInSAR) data during a single pass of the airborne platform. This advanced functionality is based on the sequential use of two sets of antenna elements having each transmit and receive capabilities. This sequence of acquisitions is illustrated in Figure 3, and it may be described as follows. A signal is first emitted by Antenna 1 at position $S_1$ and the response of the scene is received by both Antenna 1 and Antenna 2. The corresponding signals are called $y_1$ and $y_2$. Then, a signal is emitted by Antenna 2 at position $S_2$ and the scene echo is measured by the same antenna, which is called $y_3$. Let ${\varphi }_{i,j}$ represent the phase delay between the signals $y_i$ and $y_j$ due to the distance travelled by the measured wave; according to the configuration sketched in Figure 3, one can then write

$${\varphi }_{1,2}=\frac{2\pi }{\lambda }(2r_1-2r_2)=\frac{4\pi }{\lambda }(r_1-r_2),{\varphi }_{1,3}=\frac{2\pi }{\lambda }(r_1-r_2)=\frac{4\pi }{\lambda }\frac{(r_1-r_2)}{2}$$

(1)

The relationship between the components of ${\varphi }_{i,j}$ related to the acquisition geometry is implied by the specific ping-pong, non-ping-pong measurement sequence. Such an acquisition is equivalent to using a three-element equi-spaced array, with an inter-element spacing equal to $B/2$, as shown in Figure 3. The resulting dual-baseline interferometric and polarimetric SAR system may be used to reconstruct the scene 3D reflectivity, as shown in [17] for urban areas.

Compared to a classical repeat-pass acquisition system, the dual-baseline TomoSAR configuration used here presents several advantages. First, it measures a complete MB-InSAR set within a single flight, instead of three in the repeat-pass configuration, and it can significantly reduce the experimental cost and efforts. Second, single-pass acquisitions are only marginally affected by temporal decorrelation, which may reveal a very limiting factor, particularly over forested areas whose response is known to be highly affected by the presence of wind. The main drawback linked to the use of a single-pass tomographic system concerns the restricted range of baseline B, which may be physically realized on an airborne platform. In the experimental results of this paper, the baseline B was set to 3.5 m, and the corresponding vertical resolution ${\delta }_z$ varied from 13 m in near range to 84 m in far range, as shown in Figure 4.

The correlation properties of the images acquired by the system on the test site are illustrated in Figure 5. One can note that the level of correlation of the InSAR data pair corresponding to the smaller baseline, $B/2$, is characterized by high values, even over densely vegetated areas, and does not depend on the considered polarimetric channel. In the larger baseline case, B, volume decorrelation is observed in near range and becomes negligible as range increases, with a high level of correlation in mid and far range. These correlation maps are characteristic of an absence of temporal decorrelation.

3. Estimation of Forest DTM and DSM

Polarimetric tomographic processing techniques, adapted to dual-baseline interferometric and polarimetric SAR system, are applied for estimating the forest DTM, DSM and tree heights.

3.1. Polarimetric Tomographic Processing

In the case of an MB-InSAR acquisition, considering an azimuth-range resolution cell that contains $n_s$ backscattering contributions from scatterers located at different heights and assuming no decorrelation between the different acquisitions, the received data vector, $\mathsf{y}\in {\mathbb{C}}^{M\times 1}$, can be formulated as:

$$\mathsf{y}\left(l\right)=\sum _{i=1}^{n_s}{s}_i\mathsf{a}\left(z_i\right)+\mathsf{n}\left(l\right)=\mathsf{A}\left(\mathsf{z}\right)\mathsf{s}+\mathsf{n}\left(l\right)$$

(2)

where $l=1,\cdots ,L$ indicates one of the L independent realizations of the signal acquisition, also called looks. The source signal vector, $\mathsf{s}={[s_1,\cdots ,s_{n_s}]}^T$, contains the unknown complex reflection coefficient of the $n_s$ scatterers, and $\mathsf{n}\in {\mathbb{C}}^{M\times 1}$ represents the complex additive noise, assumed to be Gaussianly distributed with zero mean variance ${\sigma }_n^2$, and to be white in time and space, i.e., $\mathsf{n}\sim \mathcal{N}(0,{\sigma }_n^2{\mathsf{I}}_{(M\times M)})$ and E$\left(\mathsf{n}\left(l\right){\mathsf{n}}^{†}\left(k\right)\right)={\sigma }_n^2{\mathsf{I}}_{(M\times M)}{\delta }_{l,k}$. The steering vector $\mathsf{a}\left(z\right)$ contains the interferometric phase information associated to a source located at the elevation position z above the reference focusing plane and is given by:

$$\mathsf{a}\left(z\right)={[1,exp\left(j{k}_{z_2}z\right),\cdots ,exp\left(j{k}_{z_M}z\right)]}^T$$

(3)

where $k_{z_j}=\frac{4\pi }{\lambda }\frac{B_{{\perp }_j}}{r_{1}sin\theta }$ is the two-way vertical wavenumber between the master and the jth acquisition tracks, with the corresponding perpendicular baseline $B_{{\perp }_j}$. The carrier wavelength is represented by $\lambda$, whereas $\theta$ stands for the incidence angle and $r_1$ is the slant range distance between the master track and the scatterer. The steering matrix $\mathsf{A}\left(z\right)$ consists of $n_s$ steering vectors $\mathsf{A}\left(z\right)=[\mathsf{a}\left(z_1\right),\cdots ,\mathsf{a}\left(z_{n_s}\right)]$ with $\mathsf{z}={[z_1,\cdots ,z_{n_s}]}^T$.

The polarimetric response of a scatterer may be represented by a $N_P$-element data vector, v, containing the complex scattering coefficients measured for different transmitted and received polarization states. In a monostatic configuration, and using the well-known Pauli polarization basis [18], this vector is written as $\mathsf{v}=\frac{1}{\sqrt{2}}[S_{\mathrm{hh}}+S_{\mathrm{vv}},S_{\mathrm{hh}}-S_{\mathrm{vv}},2S_{\mathrm{hv}}]$, where $S_{\mathrm{xy}}$, with x and y equal to h or v, represents one elements of the $(2\times 2)$ scattering matrix S. Signals measured in an MB-PolInSAR configuration are then gathered in an $N_{P}M$-element polarimetric tomographic received signal, ${\mathsf{y}}_P$, formed by stacking the TomoSAR responses for each polarization channel, as ${\mathsf{y}}_P={[{\mathsf{y}}_1^T\cdots ,{\mathsf{y}}_{N_P}^T]}^T\in {\mathbb{C}}^{N_{P}M\times 1}$, where ${\mathsf{y}}_q\in {\mathbb{C}}^{N_{P}M\times 1}$ represents the MB-InSAR response for the qth polarimetric channel, i.e., ${\left[{\mathsf{y}}_q\right]}_j={\left[{\mathsf{v}}_j\right]}_q$ ($q=1,\cdots ,N_P$). Using this convention, the polarimetric steering vector of the ith source, i.e., the ideal interferometric response with unitary polarimetric target vector ${\mathsf{k}}_i$, is given by

$$\mathsf{a}(z_i,{\mathsf{k}}_i)=\left({\mathsf{I}}_{N_P\times N_P}\otimes \mathsf{a}\left(z_i\right)\right){\mathsf{k}}_i={\mathsf{B}}_{\mathsf{a}}\left(z_i\right){\mathsf{k}}_i$$

(4)

where ${\mathsf{k}}_i$ is the polarimetric target vector of the ith source and $\mathsf{K}=[{\mathsf{k}}_1,\cdots ,{\mathsf{k}}_{n_s}]$. Similar to the single polarization expression in (2), the received MB-PolInSAR data ${\mathsf{y}}_P\left(l\right)$ is given by

$${\mathsf{y}}_P\left(l\right)=\sum _{i=1}^{n_s}{s}_i\mathsf{a}(z_i,{\mathsf{k}}_i)+\mathsf{n}\left(l\right)=\mathsf{A}(\mathsf{z},\mathsf{K})\mathsf{s}\left(l\right)+\mathsf{n}\left(l\right)\in {\mathbb{C}}^{3M\times 1}$$

(5)

where $\mathsf{s}\left(l\right)\in {\mathbb{C}}^{n_s\times 1}$ represents a realization of the complex amplitude of the ith source. The reflectivity, scattering vector and height ${\sigma }_i,{\mathsf{k}}_i,z_i$ of each source, may be estimated from ${\mathsf{R}}_P=\mathsf{E}\left({\mathsf{y}}_P\left(l\right){\mathsf{y}}_P^{†}\left(l\right)\right)$, but in practice these parameters are estimated from the $L-$look estimate of the data covariance matrix ${\widehat{\mathsf{R}}}_P=\frac{1}{L}{\sum }_1^L{\mathsf{y}}_P\left(l\right){\mathsf{y}}_P^{†}\left(l\right)$ using polarimetric tomographic techniques. Compared to single polarization tomography, the use of several polarization channels results in an undeniable increase of performance improvement for recovering the physical features of the observed medium. Polarimetric diversity also provides a better estimation accuracy for polarimetric tomographic estimators than single polarization ones, as shown in [17,19]. Various polarimetric tomographic techniques are listed in [17] where their performance is compared in details.

3.2. Limitation Due to a Coarse Vertical Resolution

Natural volumes, such as forest canopies, being composed of a large number of scatterers whose responses cannot be discriminated by the resolution of analysis, their scattering patterns are generally considered as a vertical density of random or speckle-affected reflectivity. To estimate this type of power spectral density, the classical polarimetric tomographic estimator, P-Capon, is usually a good choice, and it is well appreciated for the characterization of forests [14,16], as it allows estimating the vertical reflectivity with an improved resolution with respect to the traditional Beamforming technique. It is non-parametric and estimates the reflectivity at a given elevation as:

$$\begin{array}{cc}\hfill P_C\left(z\right)& =\underset{\mathsf{k}}{max}{P}_C(z,k)=\frac{1}{{\lambda }_{min}\left({\mathsf{B}}_a^{†}\left(z\right){\widehat{\mathsf{R}}}_P^{-1}{\mathsf{B}}_{\mathsf{a}}\left(z\right)\right)}\hfill \\ \hfill {\widehat{\mathsf{k}}}_{C_i}& ={\mathsf{e}}_{min}\left({\mathsf{B}}_{\mathsf{a}}^{†}\left({\widehat{z}}_i\right){\widehat{\mathsf{R}}}_P^{-1}{\mathsf{B}}_{\mathsf{a}}\left({\widehat{z}}_i\right)\right)\hfill \end{array}$$

(6)

where ${\widehat{\mathsf{k}}}_{C_i}$ represents the polarimetric target vector at the height $z_i$, ${\lambda }_{min}$ is the smallest eigenvalue of ${\mathsf{B}}_a^{†}\left(z\right){\widehat{\mathsf{R}}}_P^{-1}{\mathsf{B}}_{\mathsf{a}}\left(z\right)$ and its corresponding eigenvector it ${\mathsf{e}}_{min}$.

As mentioned above, the dual-baseline TomoSAR configuration used in this study, possesses a vertical resolution ranging from 13 m in near range up to 84 m in far range, whereas the studied forest stands have tree heights varying from 10 to 30 m. Despite the known capability of Capon for improving the resolution of spectral analysis with respect to the classical Fourier or Beamforming method, difficulty may arise with respect to a precise characterization of the forest stands using the present TomoSAR configuration. Three test lines are respectively taken at near, mid and far range, as shown in Figure 6a, to evaluate the performance of Capon tomograms for estimating ground and tree top elevations. The resulting P-Capon tomograms in Figure 6 cannot be used to precisely characterize the forest structures in the vertical direction. Due to the resolution limitation, the profiles of vertical reflectivity over forested areas are generally composed of a single lobe merging ground and volume contributions, as shown with the normalized vertical reflectivity profiles given in Figure 7.

3.3. Proposed Solution Based on High-Resolution Spectral Analysis

High-Resolution (HR) spectral analysis techniques are able to estimate the elevation of well-localized contributions with a very high precision, that does not depend on the signal-to-noise ratio, but rather on the quality of the estimate of the data covariance matrix [20]. As illustrated in [11], HR techniques are not adapted to the characterization of distributed, continuous and speckle-affected densities of reflectivity, such as those corresponding to a forest response. Over such media, HR techniques provide estimates whose elevation is randomly distributed within the forest limits, as shown in [21]. Nevertheless, in the case of very coarse vertical resolution, it is very likely that the density of reflectivity of a forest can be approximated by two well-localized, i.e., Dirac-like, contributions. The first one is the ground response and contains the echoes directly reflected by the underlying ground as well as those corresponding to double-bounce reflection between the ground and tree trunks. The second one corresponds to the canopy, whose reflectivity profile vertical spread is much smaller than its physical extent, due to the extinction, and whose thickness remains by far inferior to the vertical resolution of the tomographic analysis. Under these assumptions, HR techniques may hence be used to estimate the ground elevation as well as the phase center of the canopy, as two well localized contributions.

Most HR polarimetric tomographic estimators are multi-dimensional approaches, which provide precise estimates of elevation and polarimetric parameters, even for coherent scatterers, as demonstrated in [17]. Nevertheless, multi-dimensional approaches usually suffer from a high computational complexity, and, to perform tomographic processing over a large forested area, a mono-dimensional HR approach is computationally attractive. The P-MUSIC technique is hence a good choice to localize the phase centers of the ground ($z_g$) and volume ($z_v$) scattering contributions in the vertical direction. Its criterion is given by:

$$\begin{array}{cc}\hfill P_M\left(z\right)& =\underset{\mathsf{k}}{max}{P}_M(z,\mathsf{k})=\frac{1}{{\lambda }_{min}({\mathsf{B}}_{\mathsf{a}}^{†}\left(z\right){\widehat{\mathsf{E}}}_n{\widehat{\mathsf{E}}}_n^{†}{\mathsf{B}}_{\mathsf{a}}\left(z\right))}\hfill \\ \hfill {\widehat{\mathsf{k}}}_{M_i}& ={\mathsf{e}}_{min}({\mathsf{B}}_{\mathsf{a}}^{†}\left({\widehat{z}}_i\right){\widehat{\mathsf{E}}}_n{\widehat{\mathsf{E}}}_n^{†}{\mathsf{B}}_{\mathsf{a}}\left({\widehat{z}}_i\right))\hfill \end{array}$$

(7)

where the noise subspace ${\mathsf{E}}_n=[{\mathsf{e}}_{n_s+1},\cdots ,{\mathsf{e}}_{N_{P}M}]$ is spanned by the eigenvectors corresponding to the $N_{P}M-n_s$ smallest eigenvalues of ${\mathsf{R}}_P$, assuming a number of sources equal to $n_s$.

To estimate the ground and volume phase centers in elevation, one may naturally set $n_s=2$. Indeed, the results given in Figure 8 indicate that, with such a parameterization, the P-MUSIC technique cannot precisely estimate the phase centers of ground and volume contributions in elevation. This is due to the fact that the covariance matrices of polarimetric responses of ground and canopy contributions have a rank larger than 1 and cannot be well described by a single target vector as in (7).

3.3.1. Model Order Selection for Ground and Volume Separation

In the general TomoSAR case, i.e., without any particular assumption of the source signal statistics, M acquisitions, composed of $N_P$ polarimetric channels, may be used to reconstruct up to $M-1$ sources having equal polarimetric responses and up to $N_{P}M-1$ sources having all different polarimetric patterns. Considering ground and canopy polarimetric covariance matrices having full rank, i.e., occupying each a subspace of dimension equal to 3, the number of selected sources should be at least 4, i.e., $n_s\ge 4$. Examples computed for $n_s>2$, in Figure 9, demonstrate that the ground and volume phase centers can be separated by the P-MUSIC estimator with $n_s=4$, in Figure 9c, whereas a better quality of tomographic focusing is provided by $n_s=5$ in Figure 9d. The model order $n_s=5$ is selected for the tomographic processing of forested areas in this paper.

To confirm this finding, a dual-baseline L-band polarimetric dataset is simulated over a homogenous forested area with varying forest parameters and ground conditions, using the PolSARproSIM software [22]. The resulting Pauli-coded images and the corresponding coherences are shown in Figure 10, with the PolSARproSIM parameters given in Figure 10a. The P-Capon technique and P-MUSIC estimator are both applied to the simulated dataset. Considering an intermediate ground roughness, i.e., Roughness = 4, the P-Capon method and P-MUSIC estimator with $n_s=2$ are able to estimate the ground and volume phase centers in the vertical direction, as shown in Figure 11b,e. Nevertheless they cannot localize the ground scattering contribution for rough surfaces (Figure 11c,d,f,g). With a model order increased to 4, the phase centers of ground and volume contributions are precisely estimated by the P-MUSIC technique, even for very rough surfaces, as illustrated in Figure 11h–j. These results from simulated data further validate the above empirical approach for model order selection.

3.3.2. Tree Height Tomographic Estimation Approach

As described above, the P-MUSIC spectral estimation technique is able to accurately estimate the ground elevation as well as the vertical position of the canopy volume phase center. To estimate the tree top elevation, a technique based on the estimate of the vertical reflectivity density provided by the P-Capon technique, formerly used in [8,14,23], and suggested by the profiles given in Figure 7, is adopted. The tree top height, $z_{top}$, is estimated as the elevation so that $P_c\left(z_{top}\right)=P_c\left(z_v\right)-3$, with $z_{top}>z_v$, as illustrated in Figure 12. The tree height $h_{tree}$ can be obtained in such a way: $h_{tree}=z_{top}-z_g$.

To demonstrate the effectiveness of the proposed approach, the ground topography and tree top heights are estimated for three test lines, as shown in Figure 13. The estimated elevations for the underlying ground and forest top show an excellent match with LiDAR-derived estimates.

This technique is then applied to the whole test site so as to map the forest DTM, DSM and tree heights.

4. Forest/Non-Forest Mapping for Model Order Selection

Forest/non-forest Mapping (FM) is generally required for GIS products related to forest characterization. The quality of this forest mapping is relevant to the final accuracy of forest-related variables. In this paper, FM is required to adequately set the model order of HR tomographic processing. Over non-forested area, the only scattering contribution to be estimated is the ground response, and ${\widehat{n}}_s=1$. Over forested areas, according to the study presented in the previous section, ${\widehat{n}}_s=5$ is used. Forest mapping is achieved using the PolInSAR correlation information between $y_1$ and $y_2$, i.e., using the larger baseline. The magnitude and phase of ${\gamma }_{{opt}_1}$, defined as the optimal PolinSAR coherence, i.e., the coherence estimated over the polarimetric channel leading to the maximum magnitude, defined by Cloude and Papathanassiou [24] and studied in detail by Ferro-Famil et al. [25], are shown in Figure 14.

The forest mask is derived by thresholding the magnitude of this optimal coherence, $|{\gamma }_{{opt}_1}|$ in the following way:

$$\mathrm{FM}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}|{\gamma }_{{opt}_1}|<T\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(8)

The threshold value $T=0.93$ is selected, and the resulting FM mask, shown in Figure 15d, is compared to the LiDAR-derived tree heights given in Figure 15b. A very good agreement is found between the FM mask of Figure 15d and the binary map of LiDAR tree heights larger than 10 m in Figure 15c. Such a map effectively corresponds to a forest/non-forest map of the test site.

5. Global Tomographic Processing over the Test Site

Based on the FM mask of Figure 15d, the model order is set to $n_s=1$ for non-forest areas and $n_s=5$ for forested areas. The proposed tomographic approach is applied over the test site of Edson. The resulting DTM and DSM are validated against LiDAR-derived measurements. The TomoSAR-estimated DTM in Figure 16a shows some boundary effects at non-forest/forest transitions, as well as a slight overestimation of the elevation of the underlying ground. Nevertheless, it globally coincides with the LiDAR-derived DTM, with a comparable high resolution. Figure 16b demonstrates the undeniable quality of estimated DSM compared to the LiDAR-derived one. Some underestimated TomoSAR-derived values indicated by red circles in Figure 16b may be due to the forest/non-forest mask errors, as suggested by Figure 16. Compared with the LiDAR-derived tree heights of Figure 16c, the proposed tomographic technique provides globally a good performance for forest height estimation, using a dual-baseline TomoSAR configuration formed by a single-pass PolInSAR system with a single physical baseline and an extremely coarse vertical resolution.

To provide a detailed quantitative evaluation of the estimated forest parameters, nine Regions of Interest (ROI) are selected over forested areas, as shown in Figure 17, each covering a 55 $\mathrm{m}$× 55 $\mathrm{m}$ wide area. Statistical features of TomoSAR- and LiDAR-derived tree heights are calculated, for each ROI, as

$$\begin{array}{cc}\hfill {\widehat{m}}_{X,i}& =\frac{1}{N_i}\sum _{q\in R_i}{\widehat{x}}_{X,q}{\widehat{\sigma }}_{X,i}=\sqrt{\frac{1}{N_i-1}{\sum }_{q\in R_i}{({\widehat{x}}_{X,q}-{\widehat{m}}_{X,i})}^2}\hfill \\ \hfill \mathrm{with}& X=Tomo,Lidar,x=h_{tree}\hfill \end{array}$$

where $N_i$ represents the number of pixels included in the ith ROI, $R_i$, with $i=1,\cdots ,9$, as listed in

Table 1. TomoSAR- and LiDAR-derived tree height statistical descriptors.

ROI	${\widehat{\mathit{m}}}_{\mathit{Tomo}}\pm {\widehat{\mathit{\sigma }}}_{\mathit{Tomo}}$ (m)	${\widehat{\mathit{m}}}_{\mathit{Lidar}}\pm {\widehat{\mathit{\sigma }}}_{\mathit{Lidar}}$ (m)
1	24.57 ± 2.59	26.39 ±0.72
2	22.40 ± 2.90	23.71 ± 0.77
3	19.32 ± 1.86	21.84 ± 0.38
4	18.10 ± 2.08	19.95 ± 0.42
5	20.58 ± 1.73	21.40 ± 0.72
6	21.68 ± 1.88	26.36 ± 0.58
7	20.80 ± 3.55	24.63 ± 0.94
8	24.33 ± 1.43	23.80 ± 1.21
9	19.54 ± 1.48	21.50 ± 0.68

. To further assess the quality of the proposed technique, the differences between TomoSAR- and LiDAR-derived estimates are characterized by their first- and second-order moments, computed as

$$\begin{array}{cc}\hfill \Delta {\widehat{x}}_q& ={\widehat{x}}_{Tomo,q}-{\widehat{x}}_{Lidar,q}\mathrm{with}x=z_g,h_{tree}\hfill \\ \hfill {\widehat{m}}_{\Delta x_i}& =\frac{1}{N_i}\sum _{q\in R_i}\Delta {\widehat{x}}_q{\widehat{\sigma }}_{\Delta x_i}=\sqrt{\frac{1}{N_i-1}{\sum }_{q\in R_i}{(\Delta {\widehat{x}}_q-{\widehat{m}}_{\Delta x_i})}^2}\hfill \end{array}$$

Figure 18a reveals the estimated DTM has a quality comparable to the LiDAR-derived one, with an average overestimation of 0.87 m for these ROIs. The estimated tree heights demonstrate a bias less than 2 m for most ROIs, except for ROI 3 and ROI 6, as depicted in Figure 18b, and still achieve an average underestimation of 2.2 m, which is less than $10\%$ of the LiDAR-derived tree heights over most ROIs. One may consider the use of a dual-baseline single-pass TomoSAR configuration, i.e., of a minimal configuration, and the accuracy of estimated forest parameters is hence well appreciated.

6. Discussion

The forest DTM, DSM and tree heights estimated by the proposed tomographic approach, as displayed in Figure 16, indicates that the FM mask plays an important role in the forest estimation process. Patches of short vegetation may be identified in the LiDAR-derived DSM in Figure 15b, but not in the TomoSAR-derived one. This is due to the coherence-derived FM mask of Figure 15d, which excludes the trees with height larger than 10 m. Boundary effects between forests and bare areas are clearly shown by the estimated DTM at the Edson area in Figure 16a. As shown in the Pauli-coded image of this test site in Figure 15a, these boundary effects mostly happen at the shadow areas where the underlying DTM cannot be reconstructed due to low SNR. It may hence be reasonable to remove the boundary areas and then interpolate DTM and DSM to improve the final estimates.

From the methodological perspective, the proposed hybrid tomographic approach provides an excellent performance for this dataset, using a model order selected in an adaptive way, as described in Section 3.3.1.

From the system perspective, unlike from other single-baseline PolInSAR systems, this specific system combines both ping-pong and non-ping-pong InSAR modes simultaneously, and it results in a dual-baseline configuration that permits tomographic processing. Moreover, this system performs the data acquisitions in a single flight, which effectively reduces the technical effort and cost for forest monitoring, and is intrinsically robust to temporal decorrelation, resulting in a better 3D reconstruction of forests.

7. Conclusions

This paper addresses forest height estimation for boreal forests at the test site of Edson in Alberta, Canada, using Intermap’s single-pass PolInSAR dataset. This particular dataset is acquired by using both ping-pong and non-ping-pong modes, which permits to form a dual-baseline TomoSAR configuration, i.e., an extreme configuration of tomographic processing. Mono-dimensional tomographic techniques are preferred in this paper, considering computational complexity of tomographic processing over a large forested area. Classical polarimetric Capon and MUSIC tomographic estimators are discussed and applied to the test data. The P-Capon technique provides a robust estimate of the vertical distribution of the back scattered power over forested areas, but it cannot recover the underlying ground elevation due to its limited resolution. The P-MUSIC technique with an appropriate model order is able to separate the volume and ground contributions in this specific low resolution configuration, but it cannot recover the tree top height. To overcome these limitations, a tomographic approach is proposed, based on P-MUSIC and P-Capon techniques, which provides excellent performance for forest parameter estimation, including DTM, DSM and tree heights. In conclusion, this paper demonstrates the undeniable potential of single-pass L-band PolInSAR system for forest height estimation, using simple tomographic techniques. Besides, this system, characterized by its single-pass configuration, can effectively fight against the temporal decorrelation of forests, leading to a potential solution for the 3D forest imaging from space.