## 1.

## Introduction

By using multiple baseline acquisitions along the direction of elevation, synthetic aperture radar tomography (TomoSAR) can extend conventional two-dimensional (2-D) SAR imaging to a three-dimensional (3-D) reconstruction. For every azimuth-range pixel, this technique can reconstruct the reflectivity function along the elevation direction (this direction is perpendicular to the azimuth and range direction) by using the formation of an additional synthetic aperture in elevation. In this way, SAR tomography can lead to a refined analysis of volume structures, namely forested areas.

Actually, the focused SAR image from multiple baseline acquisitions is the Fourier transform of the reflectivity function in the elevation direction.^{1} Therefore, SAR tomography can be considered a spectrum estimation problem. There are several spectral analysis techniques used to perform tomography, such as beamforming, Capon,^{2}^{,}^{3} multiple signal classification,^{4} truncated singular value decomposition,^{4} and weighted subspace fitting estimators.^{5} However, in practice, due to unevenly sampled acquisitions, infrequent passes over the area of interest, and limited overall acquisition baseline extent, these techniques generally bring about some quality problems in the estimation accuracy of the positions of scatterers and the resolution power such that scatterers can be discriminated. Recently, compressive sensing (CS)^{6}^{,}^{7}8.9.^{–}^{10} has been used for SAR tomography. The sparse reconstruction technique aims at using a small number of acquisitions to construct the 3-D image with a super-resolution power under the condition that the reflectivity signal is sparse along the elevation direction. The theory and practice of CS-based SAR tomography have been well developed for signals in the single range-azimuth resolution cell.^{11}12.13.^{–}^{14}

However, in the analysis of forest scenes, it is difficult to use the CS technique in the retrieval of the vertical distribution of the backscattered power, mainly because the reflectivity signal of the distributed media is rarely sparse in the object domain and the optimal sparse representations of the reflectivity signal of the distributed media must be exploited in forested areas. A good solution may be found in the work by Aguilera et al.^{15}^{,}^{16} The authors exploited suitable sparse representations of forested areas in the wavelet domain. In order to further discriminate and characterize the objects under analysis using their polarimetric responses, they extended this method to the FP case.^{17} In the process of signal recovery, the method in Ref. 17 takes advantage of the inter-signal correlations between scattering mechanisms by means of distributed compressive sensing (DCS), which enables the joint recovery of multisignal ensembles through exploiting their joint sparsity.

Furthermore, for the backscatter from the same structure, the signals from polarimetric channels also have approximately joint sparsity. Therefore, in this paper, we proposed a new wavelet-based distributed compressive sensing FP MB-InSAR tomography method (FP-WDCS TomoSAR method). In the process of signal recovery, this method mostly focuses on taking advantage of the inter-signal correlations between polarimetric channels, but not scattering mechanisms by means of DCS, because the sensors from different polarimetric channels mostly observe the same target area and the ensemble of signals can be expected to possess similar joint structures.

This paper is outlined as follows. General SAR tomographic signal models and the DCS methodology in SAR tomography are given in Sec. 2 as well as the proposed method. In Sec. 3, the datasets used are listed. In Sec. 4, the super-resolution power of our approach is analyzed and the performance of our approach is demonstrated by MB-PolInSAR datasets acquired by the ONERA SETHI airborne system over a test site in Paracou, French Guiana. Finally, conclusions are drawn in Sec. 5.

## 2.

## Methodology

## 2.1.

### Tomographic Synthetic Aperture Radar Imaging Model

Through multiple baseline synthetic aperture radar (SAR) observations over the same target area at different times and slightly different orbit positions (the elevation aperture), a stack of $N$ complex SAR datasets can be obtained. Each perpendicular baseline with respect to a master track can be calculated. Let the baselines (or the elevation aperture positions) be ${b}_{n}$. After some preparation of the stack of data (registration and phase compensation, for example), the focused complex value ${g}_{n}$ of an azimuth-range pixel $({x}_{0},{r}_{0})$ of the $n$th acquisition is

## (1)

$${g}_{n}=\underset{\mathrm{\Delta}s}{\int}\gamma (s)\text{\hspace{0.17em}}\mathrm{exp}(-j2\pi {\xi}_{n}s)\mathrm{d}s,$$^{18}${\xi}_{n}=-2{b}_{n}/(\lambda {r}_{0})$ is the spatial (elevation) frequency depending on the (more or less random) elevation aperture position ${b}_{n}$, range ${r}_{0}$, and the wavelength $\lambda $. We can see from the continuous-space system model of Eq. (1) that the multibaseline data acquisition is actually a randomly sampled Fourier transform of $\gamma (s)$. Thus, an inherent Rayleigh resolution in elevation is ${\rho}_{s}=\lambda r/(2\mathrm{\Delta}b)$, where $\mathrm{\Delta}b$ is the elevation aperture size. Through approximation by discretizing the continuous reflectivity function along $s$, in the presence of noise $\u03f5$, the discrete reflectivity function-space system model can be written as where $g$ is the measurement vector with $N$ elements ${g}_{n}$, $\mathbf{R}$ is an $N\times L$ mapping matrix with ${R}_{n\times l}=\mathrm{exp}(-j2\pi {\xi}_{n}{s}_{l})$, and $\gamma $ is the discrete reflectivity vector with $L$ elements ${\gamma}_{l}=\gamma ({s}_{l})$. ${s}_{l}$ ($l=1,\dots ,L$) denotes the discrete elevation positions. Equation (2) is a sampled discrete Fourier transform of the elevation profile $\gamma (s)$.

## 2.2.

### Covariance Matching-Based Tomographic Synthetic Aperture Radar

The distributed media, like rough ground and canopy, are characterized by a scattering response having a random behavior conferred by the speckle effect. Therefore, in forested areas, we need to use the statistical nature of the scattering response of the observed objects, which can be characterized through second-order moments of data, like the data covariance matrix, and the information of the multiple range-azimuth resolution cells must be considered.

In multiple SAR observations, we expect that the main structural characteristics of the target objects remain unchanged. This consideration can easily be approximated by tomographic acquisition within a short time. This hypothesis suggests that the temporal decorrelation of these stationary target objects is negligible. Thus, the covariance matrix of $g$ can be written out as

## (3)

$$C=E\{g{g}^{H}\}=\mathbf{R}\text{\hspace{0.17em}}\mathrm{diag}(p){\mathbf{R}}^{H}+{\sigma}^{2}I,$$^{19}Moreover, Eq. (3) can be rewritten as a linear equation: with where $1\le j,k\le N$, ${\mathit{R}}_{j}$ represents the $j$’th row of $\mathbf{R}$, and $\odot $ indicates the element-wise multiplication. In this paper, the objective of TomoSAR concerns the estimation of reflectivity profile $P$ for each azimuth-range pixel from the covariance matrix of the multibaseline measured signal $g$.

## 2.3.

### Distributed Compressive Sensing

CS builds on the ground-breaking work of Candes, Romberg, Tao, and Donoho,^{7}8.9.^{–}^{10} who showed that one can recover a signal $f\in {\mathrm{R}}^{n}$ from $m$ corrupted linear measurements $b=\mathit{A}f+z$, where $\mathit{A}$ is an $m$ by $n$ sensing matrix with $m$ typically being much smaller than $n$ and $z$ being an arbitrary and unknown vector of noise. The theory asserts that if $f$ has approximately sparse representation in one basis $\mathrm{\Psi}$, then it is indeed possible to recover $f$ from a small number of projections $b$, under the condition that $\mathit{A}$ is incoherent with $\mathrm{\Psi}$, by L1 minimization:

## (5)

$$\underset{\tilde{f}}{\mathrm{min}}\Vert \mathrm{\Psi}\tilde{f}\Vert ,\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{subject to}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\Vert \mathit{A}\tilde{f}-b\Vert}_{2}\le \u03f5,$$Since the sensors presumably observe the same target area, the ensemble of signals acquired by sensors can be expected to possess some joint structures. In this way, DCS^{20} theory extends the concept of a signal being sparse in some basis to the concept of an ensemble of signals being jointly sparse. It enables the joint recovery of multisignal ensembles through exploiting both intra- and inter-signal correlation structures. The DCS technique has been thoroughly studied and can be found in the literature.^{20}^{,}^{21} In this study, we applied the DCS method to SAR tomography. One advantage of this method is that it allows us to reduce the number of measurements needed for reconstruction.^{20} Thus, for SAR tomography, this model allows us to reduce the number of 2-D SAR images.

## 2.4.

### Fully Polarimetric Wavelet-Based Distributed Compressive Sensing Tomosynthetic Aperture Radar Method

In this section, DCS inversion techniques for multibaseline PolInSAR (MB PolInSAR) of forested areas will be introduced.

As described in Sec. 2.2, Eq. (4) expressed the tomographic sensing operation (using parallel tracks) for the signal of one azimuth-range pixel in one channel. It can be written as

## (6)

$${\mathrm{b}}_{xy}=\mathit{A}{\mathrm{p}}_{xy}+{\mathrm{z}}_{xy},xy=\mathrm{H}\mathrm{H},\mathrm{H}\mathrm{V},\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathrm{or}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{V}\mathrm{V}.$$Based on the DCS technique, we can simultaneously focus on all channels with a mixed ${\mathrm{L}}_{\mathrm{2,1}}$ minimization:^{21}

## (8)

$$\underset{\tilde{G}}{\mathrm{min}}{\Vert \mathrm{\Psi}\tilde{P}\Vert}_{2,1}\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{subject to}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\Vert \mathit{A}\tilde{P}-\tilde{B}\Vert}_{F}\le \u03f5,$$^{21}

The distribution of the backscattered power over forested areas is not sparse in the object domain^{19}^{,}^{22}23.^{–}^{24} because at least the Fourier spectrum of canopy backscatter is spread along the cross-range direction. Thus, a sparse basis needs to be found. Here, we use the wavelet basis as a sparse basis in which the profile $p$ has an approximately sparse representation.^{15}^{,}^{16}

However, for all possible practical observed forested terrains, it is very complicated to find a theory providing analysis on the performance of sparsity of the reflectivity profile in the wavelet basis. Thus, the FP-WDCS TomoSAR method may not be suitable in almost all cases. The process of the FP-WDCS TomoSAR method is shown in Fig. 1.

## 3.

## Datasets

To demonstrate the potential of the outlined approach, we used a stack of $N=6$ focused and coregistered P-band FP SAR images of the forest site of Paracou, French Guiana (Fig. 2), acquired by the ONERA SETHI airborne system on August 24, 2009, within the framework of the European Space Agency (ESA) TropiSAR 2009 campaign. The most common soils in Paracou are shallow ferralitic soils. The dominant woody species at the site include Lecythidaceae, Leguminosae, Chrysobalanaceae, and Euphorbiaceae. The tomographic data were acquired with vertical baselines of 50, 100, 150, 200, and 250 feet. The spatial resolution is approximately 1.000 m in the slant range direction and 1.245 m in the azimuth direction. The look angle varies from 25 deg to 60 deg from the near to far ranges, resulting in the Rayleigh resolution along the vertical direction varying from 21.403 to 41.523 m. Finally, all acquisitions were acquired within a 2-h period, making temporal decorrelation negligible. The flight plan is reported in Table 1 and 2. Finally , the software we used to solve the minimization problems was CVX.^{25}

## Table 1

ONERA SETHI airborne system parameters.

Parameters | |
---|---|

Wavelength | 0. 7542 m (P-Band) |

Range distance | 4905 m |

Incidence angle | 35.0614 deg |

Range resolution | 1.000 m |

Azimuth resolution | 1.245 m |

## Table 2

Flight parameters.

Number | Flight date | Baseline(m) |
---|---|---|

1 | 24/08/2009 | 0 |

2 | $-14.4879$ | |

3 | $-30.1163$ | |

4 | $-43.8343$ | |

5 | $-60.0632$ | |

6 | $-74.9683$ |

## 4.

## Numerical Experiment

## 4.1.

### Super-Resolution Power of the Compressive Sensing Framework-Based Estimator

For single-channel CS-based TomoSAR, the CS theory ensures that when the measurements are partial Fourier measurements, it is possible to recover with overwhelming probability the $K$ largest $\alpha i$ from a number of measurements satisfying the inequality:^{26}27.^{–}^{28}

Thus, for single-channel CS-based TomoSAR, the super-resolution factor ${\eta}_{\mathrm{sup}}$ is^{11}

$M$: the number of scenes

$C$: very small constant

$K$: the number of nonzero coefficients of the signal

${S}_{1}$: the illuminated scene extension in the elevation direction (range of possible elevations).

See Table 3 for details of the super-resolution power of the CS-based TomoSAR technique in our numerical experiment.

## Table 3

Super-resolution power of the compressive sensing (CS)-based synthetic aperture radar tomography (TomoSAR) technique.

Number of scatterers K | Super-resolution factor ηsup (Resolution m) (M=6, C=0.017) |
---|---|

1 | 15.260 (1.403–2.720 m) |

2 | 7.655 (2.796–5.454 m) |

3 | 5.387 (3.973–7.708 m) |

## 4.2.

### Real Data

To demonstrate the performance of the FP-WDCS TomoSAR method, this section is devoted to reporting the tomographic results of the method on the real data introduced in Sec. 3.

First, in order to demonstrate the advantages and disadvantages of the approach, we carried out tomographic processing at a fixed range distance of 550 m and for 700 contiguous azimuth positions (see the red dashed line segment in the test zone shown in Fig. 3) on the stack of six focused and coregistered P-band FP SAR images of the forest site of Paracou. The size of the estimation window used for covariance estimation was set to $39\times 39\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixels}$ (ground range, azimuth) for the method.

To facilitate the interpretation of the method’s inversion results, LIDAR-based terrain and canopy elevation estimates were used in the results. The LIDAR digital elevation model (DEM) and canopy model over Paracou were provided courtesy of CIRAD and the Guyafor project. The LIDAR measurements were resampled onto the SAR slant range azimuth coordinates and spatially averaged for application to the inversion results of the method.

For the FP-WDCS TomoSAR method, the sparsifying basis we used in this forested area was based on the Daubechies Symmlet wavelet with four vanishing moments and three levels of decomposition. As a result, we obtained slices in the azimuth and elevation directions of dimensions $871.5\times 120\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{m}}^{2}$, respectively. Figure 4 shows the three normalized tomograms and the corresponding discrete sources’ distribution from peaks of the reflectivity spectrum, which result from joint reconstruction by the FP-WDCS TomoSAR method for each polarimetric channel. Through joint recovery of multiple channels, the different elevation profiles of vegetation are present in different polarimetric channels. We can see in Fig. 4 that the HH channel mainly includes the reflectivity signal from the ground, like the one contributed by ground-trunk scatterers, due to sensitivity to double-bounce reflections over dihedral-like objects. Furthermore, the reflectivity signal from the canopy is mainly present in the HV channel because it is more sensitive to volumetric scatterers.

Moreover, the elevation profiles in the HV and VV channels (see Fig. 4) show that the spatial structures of the SM associated with the canopy backscatter agree with LIDAR measurements (see the top black line in Fig. 4). However, there are obvious errors between the spatial structures of the SM from ground scattering and LIDAR measurements (see the bottom black line in Fig. 4). One possible reason is that the SM associated with ground scattering is probably dominated by trunk-ground scattering. Thus, the reflectivity signal from ground scattering is an approximately point-like reflection which has no sparse representation in a wavelet basis.

To further demonstrate the performance of the FP-WDCS TomoSAR method, two traditional nonparametric estimators, namely, the FP beamforming method and FP Capon method, were used to compare to it. The specific processes of the FP Capon method and FP beamforming method are introduced in Refs. 7,2930.–31.

Figure 5 shows the estimated power with the FP Capon and FP beamforming methods in the azimuth-elevation plane. In Fig. 5, the FP Capon estimator has a higher resolution than FP beamforming. However, the resolution limitation of both methods blurs the power profile for each pixel, meaning that the limited spectral resolution of both methods affects the accuracy of the height estimation.

Moreover, we calculated the root-mean-square error (RMSE) of ground, canopy, and tree height estimation based on LIDAR measurements with the three tomography methods (Table 4). Table 4 shows that the FP-WDCS TomoSAR method has the best performance in the RMSE of ground, canopy, and tree height estimation compared with the other two methods.

## Table 4

Root-mean-square error (RMSE) of height estimation with four TomoSAR methods. Mean of tree height is 28.4790 m.

RMSE | FP-WDCS | FP-beamforming | FP-Capon |
---|---|---|---|

Ground height | 1.7182 | 3.2284 | 2.3563 |

Canopy height | 3.6669 | 4.8317 | 5.0071 |

Tree height | 4.5389 | 6.3499 | 5.9288 |

## 5.

## Conclusions and Future Work

In this paper, we proposed a new DCS framework-based FP MB-InSAR tomography method that allows a jointly reflectivity profile reconstruction from FP channels in forested areas along the elevation direction using signals acquired from a small number of orbits and different polarization channels.

Based on the assumption that the reflectivity signal of forested areas along the elevation direction is sparse or compressible in a wavelet domain, the FP-WDCS TomoSAR method enables the joint recovery of FP signal ensembles through the DCS technique to exploit both intra- and inter-signal correlation structures between polarimetric channel signals.

SP analysis shows clearly that the CS framework-based estimators provide a better super-resolution power for localizing objects with complex structures compared to features with a limited resolution shown by the traditional nonparametric estimators. Real data experiments show that the estimation accuracy of ground and canopy height obtained using the proposed method is much better than the traditional nonparametric estimators.

Therefore, further work will focus on trying to compare and analyze the performance of the estimation accuracy of canopy height using the DCS method with or without sparse expansions in the wavelet domain. Moreover, a detailed evaluation of the estimation accuracy and super-resolution power and robustness of those techniques through simulation experiments will be carried out.

## Acknowledgments

The authors gratefully acknowledge the support for this research from the Major International Cooperation and Exchange Project of National Natural Science Foundation of China (Grant No. 41120114001), and the director innovation foundation project of CEODE, CAS (No. Y2ZZ20101B). The authors would like to thank ESA EO Campaigns data Project for providing the SAR tomographic dataset. The authors also would like to thank Principal Investigator Lilian Blanc, Mr. Bruno Herault, and Mr. Grégoire Vincent with CIRAD for kindly providing the LiDAR DEMs and canopy model over Paracou.

## References

## Biography

**Lei Liang** is currently working toward his PhD in microwave remote sensing at the Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing, China. He is also with the University of Chinese Academy of Sciences, Beijing. He received his MSc degree in computational mathematics from Beijing Jiaotong University, Beijing, China, in 2010. His research interests include polar remote sensing, microwave radiometer remote sensing, synthetic aperture radar interferometry, SAR polarimetry, and SAR tomography.

**Xinwu Li** is a professor in the Institute of Remote Sensing and Digital Earth (RADI), Chinese Academy of Sciences (CAS). He received his PhD from the Institute of Remote Sensing Applications, Chinese Academy of Sciences (CAS) in 2002. He is the author of more than 17 SCI journal papers. His current research interests include microwave remote sensing, global environmental change, urban environmental remote sensing, and so on.

**Xizhang Gao** is an assistant professor in the Institute of Geographical Sciences and Natural Resources Research, Chinese Academy of Sciences. He received his PhD from Zhejiang University in 2004. His current research interests include new techniques and applications for geography information systems, spatial databases, land use classification, eco-environment quality evaluation, environmental remote sensing, and so on.

**Huadong Guo** is a full-time professor with the Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing, China. He received his MSc degree in remote sensing from the Institute of Remote Sensing Application, Chinese Academy of Sciences, Beijing, China, in 1981.Currently, he is in charge of three projects related to microwave remote sensing and space technology for global environmental change. His research interests include radar remote sensing information mechanism and applications.