<figure><img src="image_1.png" /></figure>

HAL
open science

# Forest height estimation using a single-pass airborne L-band polarimetric and interferometric SAR system and tomographic techniques

Yue Huang, Q. Zhang, Laurent Ferro-Famil

► To cite this version:

Yue Huang, Q. Zhang, Laurent Ferro-Famil. Forest height estimation using a single-pass airborne L-band polarimetric and interferometric SAR system and tomographic techniques. Remote Sensing, 2021, 13 (3), 10.3390/rs13030487. [https://hal.science/hal-03158184](https://hal.science/hal-03158184)

HAL Id: hal-03158184

https://hal.science/hal-03158184

Submitted on 25 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
<figure><img src="image_2.png" /></figure>

Article

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

Yue Huang <fnref n="1" />*, Qiaoping Zhang <fnref n="2" /> and Laurent Ferro-Famil <fnref n="1" />,<fnref n="3" />

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

**Keywords:** forest height estimation; SAR Tomography; single-pass InSAR

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Huang, Y.; Zhang, Q.; Ferro-Famil, L. Forest Height Estimation Using a Single-Pass Airborne L-Band Polarimetric and Interferometric SAR System and Tomographic Techniques. *Remote Sens.* **2021**, *13*, 487.
https://doi.org/10.3390/rs13030487

Received: 31 December 2020

Accepted: 26 January 2021

Published: 30 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

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

<footnote n="1">¹ IETR Laboratory, F-35000 Rennes, France; laurent.ferro-famil@univ-rennes1.fr</footnote>

<footnote n="2">² Intermap Technologies Corporation, Calgary, AB T3E 6K6, Canada; qzhang@intermap.com</footnote>

<footnote n="3">³ CESBIO Laboratory, 31065 Toulouse, France</footnote>

<footnote>* Correspondence: yue.huang@univ-rennes1.fr</footnote>
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].
<figure><img src="image_5.png" /><figcaption>(a) Geographic location indicated by a red box in Google Earth image (courtesy Google Earth)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Airphoto (courtesy Valtus Imagery) surrounded by ground photos</figcaption></figure>

<caption>**Figure 1.** Illustration of the Edson test site.</caption>
<figure><img src="image_7.png" /></figure>

<table><caption>**Figure 2.** LiDAR-derived tree heights over Edson (courtesy Terrapoint Canada).</caption><thead><tr><th>Region</th><th>Type</th><th>Range</th><th>LiDAR DSM-bald DTM<br>Mean ± Std (m)</th></tr></thead><tbody><tr><td>F1</td><td>Forest</td><td>Mid-Far</td><td>22.82 ± 1.21</td></tr><tr><td>F2</td><td>Forest</td><td>Mid-Far</td><td>24.59 ± 1.51</td></tr><tr><td>F3</td><td>Forest</td><td>Far</td><td>23.27 ± 1.41</td></tr><tr><td>F4</td><td>Forest</td><td>Far</td><td>22.39 ± 0.81</td></tr><tr><td>F5</td><td>Forest</td><td>Far</td><td>23.32 ± 1.26</td></tr><tr><td>F6</td><td>Forest</td><td>Mid-Far</td><td>27.74 ± 1.01</td></tr><tr><td>F7</td><td>Forest</td><td>Far</td><td>26.52 ± 0.96</td></tr><tr><td>F8</td><td>Forest</td><td>Far</td><td>28.77 ± 0.80</td></tr></tbody></table>

## 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 Intersmap'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 $\phi_{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

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

The relationship between the components of $\phi_{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.

<figure><img src="image_8.png" /><figcaption>Figure 3. Illustration of transmit-receive modes and equivalent dual-baseline TomoSAR configuration. $S_1$ and $S_2$ are the real positions of Antenna 1 and Antenna 2, whereas $S_3$ is the virtual resulting position.</figcaption></figure>

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.

<figure><img src="image_9.png" /><figcaption>Figure 4. Tomographic resolution $\delta_z$ w.r.t range.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 5. InSAR coherences over the test site.</figcaption></figure>
### 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, $\mathbf{y} \in \mathbb{C}^{M \times 1}$, can be formulated as:

$$ \mathbf{y}(l) = \sum_{i=1}^{n_s} s_i \mathbf{a}(z_i) + \mathbf{n}(l) = \mathbf{A}(\mathbf{z})\mathbf{s} + \mathbf{n}(l) \quad (2) $$

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

$$ \mathbf{a}(z) = [\exp(jk_{z_2}z), \dots, \exp(jk_{z_M}z)]^T \quad (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 $\mathbf{A}(\mathbf{z})$ consists of $n_s$ steering vectors $\mathbf{A}(\mathbf{z}) = [\mathbf{a}(z_1), \dots, \mathbf{a}(z_{n_s})]$ with $\mathbf{z} = [z_1, \dots, z_{n_s}]^T$.

The polarimetric response of a scatterer may be represented by a $N_p$-element data vector, $\mathbf{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 $\mathbf{v} = \frac{1}{\sqrt{2}}[S_{hh} + S_{vv}, S_{hh} - S_{vv}, 2S_{hv}]$, where $S_{xy}$, with x and y equal to h or v, represents one elements of the $(2 \times 2)$ scattering matrix $\mathbf{S}$. Signals measured in an MB-PolInSAR configuration are then gathered in an $N_p M$-element polarimetric tomographic received signal, $\mathbf{y}_p$, formed by stacking the TomoSAR responses for each polarization channel, as $\mathbf{y}_p = [\mathbf{y}_1^T \dots, \mathbf{y}_{N_p}^T]^T \in \mathbb{C}^{N_p M \times 1}$, where $\mathbf{y}_q \in \mathbb{C}^{N_p M \times 1}$ represents the MB-InSAR response for the qth polarimetric channel, i.e., $[\mathbf{y}_q]_q = [\mathbf{v}_j]_q$ ($q = 1, \dots, N_p$). Using this convention, the polarimetric steering vector of the ith source, i.e., the ideal interferometric response with unitary polarimetric target vector $\mathbf{k}_i$, is given by

$$ \mathbf{a}(z_i, \mathbf{k}_i) = (\mathbf{I}_{N_p \times N_p} \otimes \mathbf{a}(z_i)) \mathbf{k}_i = \mathbf{B}_a(z_i) \mathbf{k}_i \quad (4) $$

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

$$ \mathbf{y}_p(l) = \sum_{i=1}^{n_s} s_i \mathbf{a}(z_i, \mathbf{k}_i) + \mathbf{n}(l) = \mathbf{A}(\mathbf{z}, \mathbf{K})\mathbf{s}(l) + \mathbf{n}(l) \in \mathbb{C}^{3M \times 1} \quad (5) $$

where $\mathbf{s}(l) \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$, $\mathbf{k}_i$, $z_i$ of each source, may be estimated from $\mathbf{R}_p = E(\mathbf{y}_p(l)\mathbf{y}_p^\dagger(l))$, but in practice these parameters are estimated from the L-look estimate of the data covariance matrix $\hat{\mathbf{R}}_p = \frac{1}{L} \sum_{l=1}^L \mathbf{y}_p(l)\mathbf{y}_p^\dagger(l)$ 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{aligned}
P_C(z) &= \max_{\mathbf{k}} P_C(z, \mathbf{k}) = \frac{1}{\lambda_{\min}(\mathbf{B}_a^\dagger(z) \hat{\mathbf{R}}_p^{-1} \mathbf{B}_a(z))} \\
\hat{\mathbf{k}}_{C_i} &= \mathbf{e}_{\min}(\mathbf{B}_a^\dagger(\hat{z}_i) \hat{\mathbf{R}}_p^{-1} \mathbf{B}_a(\hat{z}_i))
\end{aligned}
\qquad (6) $$

where $\hat{\mathbf{k}}_{C_i}$ represents the polarimetric target vector at the height $z_i$, $\lambda_{\min}$ is the smallest eigenvalue of $\mathbf{B}_a^\dagger(z) \hat{\mathbf{R}}_p^{-1} \mathbf{B}_a(z)$ and its corresponding eigenvector it $\mathbf{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.
<figure><img src="image_11.png" /><figcaption>(a) Test lines</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) L1 ($\delta_z$ = 25.2 m)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) L2 ($\delta_z$ = 42.3 m)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) L3 ($\delta_z$ = 66.0 m)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 6.** P-Capon tomograms for test lines located in near, mid and far range, respectively. The black line is the LiDAR ground and tree top elevation estimates.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(a) Around Pixel P₁ of L1</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Around Pixel P₂ of L2</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) Around Pixel P₃ of L3</figcaption></figure>

<caption>**Figure 7.** Normalized vertical reflectivity profiles estimated by P-Capon technique at the locations indicated in Figure 6. Elevation here is measured with respect to the LiDAR ground elevation.</caption>

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{aligned}
P_M(z) &= \max_{\mathbf{k}} P_M(z, \mathbf{k}) = \frac{1}{\lambda_{\min}(\mathbf{B}_a^\dagger(z) \hat{\mathbf{E}}_n \hat{\mathbf{E}}_n^\dagger \mathbf{B}_a(z))} \\
\hat{\mathbf{k}}_{M_i} &= \mathbf{e}_{\min}(\mathbf{B}_a^\dagger(\hat{z}_i) \hat{\mathbf{E}}_n \hat{\mathbf{E}}_n^\dagger \mathbf{B}_a(\hat{z}_i))
\end{aligned} 
\quad (7) $$

where the noise subspace $\mathbf{E}_n = [\mathbf{e}_{n+1}, \dots, \mathbf{e}_{N_p M}]$ is spanned by the eigenvectors corresponding to the $N_p M - n_s$ smallest eigenvalues of $\mathbf{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).

<figure><img src="image_19.png" /><figcaption>Figure 8. P-MUSIC tomograms for test lines located in near, mid and far range, respectively. $n_s = 2$. The black line is the LiDAR ground and tree top elevation estimates.</figcaption></figure>

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

<figure><img src="image_20.png" /><figcaption>**Figure 9.** P-MUSIC tomograms with model order larger than 2. The black line is the LiDAR ground and tree top elevation estimates.</figcaption></figure>

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,i,j. These results from simulated data further validate the above empirical approach for model order selection.
<figure><img src="image_21.png" /><figcaption>**Figure 10.** Simulated dual-baseline FP data of forests. The illustrated Pauli-coded images are generated in the case of Roughness = 4.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 11.** P-Capon and P-MUSIC tomograms with $n_s = 2$ and $n_s = 4$ over a test line for different ground roughness: P-Capon (b–d); and P-MUSIC (e–j).</figcaption></figure>

### 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(z_{top}) = P_c(z_v) - 3$ dB, 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_v$.
<figure><img src="image_23.png" /><figcaption>**Figure 12.** Proposed tomographic approaches for estimation tree height.</figcaption></figure>

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.

<figure><img src="image_24.png" /><figcaption>**Figure 13.** Ground topography and tree top heights estimated by the proposed tomographic approach. Background, P-Capon tomograms; gray line, estimated ground and tree top elevations; black line, LiDAR ground and tree top elevation estimates.</figcaption></figure>

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 $\hat{n}_s = 1$. Over forested areas, according to the study presented in the previous section, $\hat{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.

<figure><img src="image_25.png" /><figcaption>Figure 14. $\gamma_{opt_1}$ map over the test site.</figcaption></figure>

The forest mask is derived by thresholding the magnitude of this optimal coherence, $|\gamma_{opt_1}|$ in the following way:

$$FM = \begin{cases} 1 & \text{if } |\gamma_{opt_1}| < T \\ 0 & \text{otherwise} \end{cases} \qquad (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.

<figure><img src="image_26.png" /><figcaption>Figure 15. Generated forest/non-forest mask.</figcaption></figure>

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 m × 55 m wide area. Statistical features of TomoSAR- and LiDAR-derived tree heights are calculated, for each ROI, as

$$ \hat{m}_{X,i} = \frac{1}{N_i} \sum_{q \in R_i} \hat{x}_{X,q} \quad \hat{\sigma}_{X,i} = \sqrt{\frac{1}{N_i - 1} \sum_{q \in R_i} (\hat{x}_{X,q} - \hat{m}_{X,i})^2} \\ \text{with } X = \text{Tomo, Lidar, } x = h_{\text{tree}} $$

where $N_i$ represents the number of pixels included in the ith ROI, $R_i$, with $i = 1, ..., 9$, as listed in Table 1. 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

$$ \Delta \hat{x}_q = \hat{x}_{\text{Tomo},q} - \hat{x}_{\text{Lidar},q} \quad \text{with } x = z_g, h_{\text{tree}} \\ \hat{m}_{\Delta x_i} = \frac{1}{N_i} \sum_{q \in R_i} \Delta \hat{x}_q \quad \hat{\sigma}_{\Delta x_i} = \sqrt{\frac{1}{N_i - 1} \sum_{q \in R_i} (\Delta \hat{x}_q - \hat{m}_{\Delta x_i})^2} $$

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.
<figure><img src="image_27.png" /><figcaption>**Figure 16.** DSM, DTM and tree height maps estimated by the proposed tomographic technique and using LiDAR data.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Figure 17. Location of ROIs.</figcaption></figure>

<table><caption>Table 1. TomoSAR- and LiDAR-derived tree height statistical descriptors.</caption><thead><tr><th>ROI</th><th>m̂<sub>Tomo</sub> ± σ̂<sub>Tomo</sub> (m)</th><th>m̂<sub>Lidar</sub> ± σ̂<sub>Lidar</sub> (m)</th></tr></thead><tbody><tr><td>1</td><td>24.57 ± 2.59</td><td>26.39 ± 0.72</td></tr><tr><td>2</td><td>22.40 ± 2.90</td><td>23.71 ± 0.77</td></tr><tr><td>3</td><td>19.32 ± 1.86</td><td>21.84 ± 0.38</td></tr><tr><td>4</td><td>18.10 ± 2.08</td><td>19.95 ± 0.42</td></tr><tr><td>5</td><td>20.58 ± 1.73</td><td>21.40 ± 0.72</td></tr><tr><td>6</td><td>21.68 ± 1.88</td><td>26.36 ± 0.58</td></tr><tr><td>7</td><td>20.80 ± 3.55</td><td>24.63 ± 0.94</td></tr><tr><td>8</td><td>24.33 ± 1.43</td><td>23.80 ± 1.21</td></tr><tr><td>9</td><td>19.54 ± 1.48</td><td>21.50 ± 0.68</td></tr></tbody></table>

<figure><img src="image_29.png" /><figcaption>Figure 18. Statistical features of TomoSAR- and LiDAR-derived elevation estimates over different ROIs. The red bar has a length equal to $2\hat{\sigma}_{\Delta x_i}$, the black cross mark indicates $\hat{m}_{\Delta x_i}$, and the green bar has a length equal to $2\hat{\sigma}_{Lidar,i}$, with $i = 1, \dots, 9$.</figcaption></figure>

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

**Author Contributions:** Conceptualization, methodology and validation, Y.H.; preparation of SAR data and LiDAR measurements, validation and review, Q.Z.; and scientific discussion, review and editing, L.F.-F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded in part by the Alberta Innovates Technologies Future-R&D Associates program and in part by ANR (the French National Research Agency) and DGA (Direction Générale de l'Armement) under the ALYS project ANR-15-ASTR-0002.

**Acknowledgments:** The authors would like to thank Terrapoint Canada Inc. for providing the LiDAR ground truth data, and Valtus Imagery Services for providing the air photos.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## Abbreviations

The following abbreviations are used in this manuscript:

PolInSAR Polarimetric Interferometric SAR

TomoSAR SAR Tomography

DTM Digital Terrain Model

DSM Digital Surface Model
## References

1. Frolking, S.; Palace, M.W.; Clark, D.B.; Chambers, J.Q.; Shugart, H.H.; Hurtt, G.C. Forest disturbance and recovery: A general review in the context of spaceborne remote sensing of impacts on aboveground biomass and canopy structure. *J. Geophys. Res. Biogeosci.* **2009**, *114*. doi:10.1029/2008JG000911.

2. Treuhaft, R.; Siqueira, P. Vertical structure of vegetated land surfaces from interferometric and polarimetric radar. *Radio Sci.* **2000**, *35*, 141–178.

3. Papathanassiou, K.P.; Cloude, S.R. Single baseline polarimetric SAR interferometry. *IEEE Trans. Geosci. Remote Sens.* **2001**, *39*, 2352–2363.

4. Cloude, S.R.; Papathanassiou, K.P. Three-stage inversion process for polarimetric SAR interferometry. *IEE Proc. Radar Sonar Navig.* **2003**, *150*, 125–134.

5. Tewaldini, S. Forest Structure Retrieval from Multi-Baseline SARs. In *Remote Sensing of Biomass*; Fatoyinbo, T., Ed.; IntechOpen: Rijeka, Croatia, 2012; Chapter 2.

6. Pardini, M.; Papathanassiou, K. On the Estimation of Ground and Volume Polarimetric Covariances in Forest Scenarios with SAR Tomography. *IEEE Geosci. Remote Sens. Lett.* **2017**, *14*, 1860–1864.

7. Ferro-Famil, L.; Huang, Y.; Pottier, E. Principles and Applications of Polarimetric SAR Tomography for the Characterization of Complex Environments. In *VIII Hotine-Marussi Symposium on Mathematical Geodesy*; Springer: Berlin/Heidelberg, Germany, 2015; pp. 1–13. doi:10.1007/1345_2015_12.

8. Huang, Y.; Ferro-Famil, L.; Lardeux, C. Polarimetric SAR tomography of tropical forests at P and L-band. In Proceedings of the 2011 IEEE International Geoscience and Remote Sensing Symposium, Vancouver, BC, Canada, 24–29 July 2011; pp. 1373–1376.

9. Li, X.; Liang, L.; Guo, H.; Huang, Y. Compressive Sensing for Multibaseline Polarimetric SAR Tomography of Forested Areas. *IEEE Trans. Geosci. Remote Sens.* **2016**, *54*, 153–166.

10. Aghababaei, H.; Ferraioli, G.; Ferro-Famil, L.; Huang, Y.; Mariotti D'Alessandro, M.; Pascazio, V.; Schirinzi, G.; Tewaldini, S. Forest SAR Tomography: Principles and Applications. *IEEE Geosci. Remote Sens. Mag.* **2020**, *8*, 30–45.

11. Huang, Y.; Levy-Vehel, J.; Ferro-Famil, L.; Reigber, A. Three-Dimensional Imaging of Objects Concealed Below a Forest Canopy Using SAR Tomography at L-Band and Wavelet-Based Sparse Estimation. *IEEE Geosci. Remote Sens. Lett.* **2017**, *14*, 1454–1458.

12. Lee, S.; Pardini, M.; Kugler, F.; Papathanassiou, K.; Hajnsek, I. PolInSAR forest height inversion by means of L-band F-SAR data. *PolInSAR Workshop 2013*.

13. Hajnsek, I.; Kugler, F.; Lee, S.K.; Papathanassiou, K.P. Tropical-Forest-Parameter Estimation by Means of Pol-InSAR: The INDREX-II Campaign. *IEEE Trans. Geosc. Remote Sens.* **2009**, *47*, 481–493.

14. Tewaldini, S.; Rocca, F. Multibaseline Polarimetric SAR Tomography of a Boreal Forest at P- and L-Bands. *IEEE Trans. Geosci. Remote Sens.* **2012**, *50*, 232–246. doi:10.1109/TGRS.2011.2159614.

15. Mercer, B.; Zhang, Q.; Schwaebisch, M.; Denbina, M. 3D Topography and Forest Recovery from An L-Band Single-Pass Airborne PolInSAR System. In Proceedings of the 2009 IEEE International Geoscience and Remote Sensing Symposium, Cape Town, South Africa, 12–17 July 2009.

16. Tewaldini, S. Algebraic Synthesis of Forest Scenarios From Multibaseline PolInSAR Data. *IEEE Trans. Geosc. Remote Sens.* **2009**, *47*, 4132–4142.

17. Huang, Y.; Ferro-Famil, L. 3-D Characterization of Urban Areas Using High-Resolution Polarimetric SAR Tomographic Techniques and a Minimal Number of Acquisitions. *IEEE Trans. Geosci. Remote Sens.* **2020**, 1–18. doi:10.1109/TGRS.2020.3044150.

18. Lee, J.S.; Pottier, E. *Polarimetric Radar Imaging: From Basics to Applications*; CRC Press: Chicago, IL, USA, 2008.

19. Huang, Y.; Ferro-Famil, L.; Reigber, A. Under-Foliage Object Imaging Using SAR Tomography and Polarimetric Spectral Estimators. *IEEE Trans. Geosci. Remote Sens.* **2012**, *50*, 2213–2225. doi:10.1109/TGRS.2011.2171494.

20. Stoica, P.; Nehorai, A. MUSIC, Maximum likelihood, and Cramer-Rao Bound: Further Results and Comparisons. *IEEE Trans. Acoust. Speech Signal Process.* **1990**, *38*, 2140–2150.

21. Nannini, M.; Scheiber, R.; Moreira, A. Estimation of the Minimum Number of Tracks for SAR Tomography. *IEEE Trans. Geosc. Remote Sens.* **2009**, *47*, 531–543.

22. Williams, M.; Pottier, E.; Ferro-Famil, L.; Allain, S.; Cloude, S.; Hajnsek, I.; Papathanassiou, K.; Moreira, A.; Minchella, A.; Desnos, Y.L. Forest Coherent SAR Simulation within PolSARPro: An Educational Toolbox for PolSAR and PolInSAR Data Processing. In Proceedings of the Asian Conference on Remote Sensing, Kuala Lumpur, Malaysia, 12–16 November 2007.

23. Zhang, Q.; Huang, Y.; Schwaebisch, M.; Mercer, B.; Wei, M. Forest height estimation using single-pass dual-baseline L-band PolInSAR Data. In Proceedings of the 2012 IEEE international geoscience and remote sensing symposium, Munich, Germany, 22–27 July 2012.

24. Cloude, S.R.; Papathanassiou, K.P. Polarimetric SAR interferometry. *IEEE Trans. Geosci. Remote Sens.* **1998**, *36*, 1551–1565.

25. Ferro-Famil, L.; Lee, J.; Pottier, E. Unsupervised Classification of Natural Scenes from Polarimetric Interferometric SAR Data. In *Frontiers of Remote Sensing Information Processing*; World Scientific Publishing: Singapore, 2003.