Multibaseline polarimetric synthetic aperture radar tomography of forested areas using wavelet-based distribution compressive sensing

Abstract. The three-dimensional (3-D) structure of forests, especially the vertical structure, is an important parameter of forest ecosystem modeling for monitoring ecological change. Synthetic aperture radar tomography (TomoSAR) provides scene reflectivity estimation of vegetation along elevation coordinates. Due to the advantages of super-resolution imaging and a small number of measurements, distribution compressive sensing (DCS) inversion techniques for polarimetric SAR tomography were successfully developed and applied. This paper addresses the 3-D imaging of forested areas based on the framework of DCS using fully polarimetric (FP) multibaseline SAR interferometric (MB-InSAR) tomography at the P-band. A new DCS-based FP TomoSAR method is proposed: a new wavelet-based distributed compressive sensing FP TomoSAR method (FP-WDCS TomoSAR method). The method takes advantage of the joint sparsity between polarimetric channel signals in the wavelet domain to jointly inverse the reflectivity profiles in each channel. The method not only allows high accuracy and super-resolution imaging with a low number of acquisitions, but can also obtain the polarization information of the vertical structure of forested areas. The effectiveness of the techniques for polarimetric SAR tomography is demonstrated using FP P-band airborne datasets acquired by the ONERA SETHI airborne system over a test site in Paracou, French Guiana.


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. 1Therefore, 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. 5However, 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.][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. 17In 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 waveletbased 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.

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 nth acquisition is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 1 8 9 where γðsÞ represents the reflectivity function along elevations s and Δs describes the range of possible elevations. 18ξ n ¼ −2b n ∕ðλ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 λ.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 γðsÞ.Thus, an inherent Rayleigh resolution in elevation is ρ s ¼ λr∕ð2ΔbÞ, where Δb is the elevation aperture size.Through approximation by discretizing the continuous reflectivity function along s, in the presence of noise ε, the discrete reflectivity function-space system model can be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 7 1 1 where g is the measurement vector with N elements g n , R is an N × L mapping matrix with R n×l ¼ expð−j2πξ n s l Þ, and γ is the discrete reflectivity vector with L elements γ l ¼ γðs l Þ. s l (l ¼ 1; : : : ; L) denotes the discrete elevation positions.Equation ( 2) is a sampled discrete Fourier transform of the elevation profile γðsÞ.

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 where Ef•g is the expectation operator, (•) denotes the conjugate transpose, p is a non-negative real vector, diagðpÞ is a matrix whose main diagonal equals p and is made up of zeros in its offdiagonal entries, and σ 2 is the unknown noise power. 19Moreover, Eq. ( 3) can be rewritten as a linear equation: ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 3 9 6 with E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 2 ; 1 1 6 ; 3 5 3 ; t e m p : i n t r a l i n k -; s e c 2 . 2 ; 1 1 6 ; 3 1 4 ; t e m p : i n t r a l i n k -; s e c 2 . 2 ; 1 1 6 ; 2 7 4 where 1 ≤ j; k ≤ N, R j represents the j'th row of R, and ⊙ 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.

Distributed Compressive Sensing
CS builds on the ground-breaking work of Candes, Romberg, Tao, and Donoho, 7-10 who showed that one can recover a signal f ∈ R n from m corrupted linear measurements b ¼ Af þ z, where 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 Ψ, then it is indeed possible to recover f from a small number of projections b, under the condition that A is incoherent with Ψ, by L1 minimization: where ε is an upper bound on the noise level.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,21In 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. 20Thus, for SAR tomography, this model allows us to reduce the number of 2-D SAR images.

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 In our method, we suppose that the signals throughout polarimetric channels share, approximately, the same joint structure, because we are expecting backscatter from the same forested scene.In this way, considering an ensemble of signals of all channels, it will be jointly expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 4 0 7 where ½b HH b HV b VV ¼ B, ½p HH p HV p VV ¼ P.
Based on the DCS technique, we can simultaneously focus on all channels with a mixed L 2;1 minimization: 21 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 7 3 5 min G kΨ Pk 2;1 subject to kA P − Bk F ≤ ε; where Ψ is a sparse basis and k • k 2;1 is the mixed norm (sum of the L 2 norms of the rows of a matrix).The L 2;1 norm not only promotes sparsity along rows and minimizes the energy along columns, but also exponentially reduces the probability of recovery failure in the number of columns of P. 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,16owever, 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.

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.  1 and 2. Finally , the software we used to solve the minimization problems was CVX.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 αi from a number of measurements satisfying the inequality: [26][27][28] E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 4 . 1 ; 1 1 6 ; 3 0 7 M ≥ CK log ðNÞ 4 : Thus, for single-channel CS-based TomoSAR, the super-resolution factor η sup is 11 See Table 3 for details of the super-resolution power of the CS-based TomoSAR technique in our numerical experiment.

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 × 39 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 × 120 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,29-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.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 intraand inter-signal correlation structures between polarimetric channel signals.
SP analysis shows clearly that the CS framework-based estimators provide a better superresolution 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 superresolution power and robustness of those techniques through simulation experiments will be carried out.
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.
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
Super-Resolution Power of the Compressive Sensing Framework-Based Estimator

e m p : i n t r a l i n k -; s e c 4
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).

Fig. 4
Fig. 4 (a, c, e) Reflectivity tomograms in the HH channel (a, b), (c, d), HV channel, and (e, f) VV channel.The panels have been normalized, such that the sum the height is unitary.The lower and upper black lines, respectively, represent the position of ground and top of tree in the elevation direction provided by LiDAR measurements, and (b, d, f) the lower and upper red lines are the corresponding discrete source distributions from ground and top of tree.

Fig. 5
Fig. 5 (a, c) Reflectivity tomograms of (a, b) the FP beamforming and (c, d) FP Capon methods.The panels have been normalized, such that the sum along the height is unitary.The lower and upper black lines, respectively, represent the position of ground and top of tree in the elevation direction provided by LiDAR measurements, and (b, d) the lower and upper red lines are the corresponding discrete source distributions from ground and top of tree.

Table 2
Flight parameters.

Table 1
ONERA SETHI airborne system parameters.

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

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