EVALUATION OF AIRBORNE L-BAND MULTI-BASELINE POL-INSAR FOR DEM EXTRACTION BENEATH FOREST CANOPY

DEM beneath forest canopy is difficult to extract with optical stereo pairs, InSAR and Pol-InSAR techniques. Tomographic SAR (TomoSAR) based on different penetration and view angles could reflect vertical structure and ground structure. This paper aims at evaluating the possibility of TomoSAR for underlying DEM extraction. Airborne L-band repeat-pass Pol-InSAR collected in BioSAR 2008 campaign was applied to reconstruct the 3D structure of forest. And sum of kronecker product and algebraic synthesis algorithm were used to extract ground structure, and phase linking algorithm was applied to estimate ground phase. Then Goldstein cut-branch approach was used to unwrap the phases and then estimated underlying DEM. The average difference between the extracted underlying DEM and Lidar DEM is about 3.39 m in our test site. And the result indicates that it is possible for underlying DEM estimation with airborne L-band repeat-pass TomoSAR technique.


INTRODUCTION
Digital elevation model (DEM) extraction in areas without forest or with sparse trees is the advantage of optical stereo pairs (Alrousan et al. 1997, D'Angelo et al. 2008), InSAR (Zhang and Yue 2017, Geymen 2014, Gao et al. 2017)and Pol-InSAR (Fu et al. 2016, Fu et al. 2017, Shimada, Natsuaki and Hirose 2018).While, DEM beneath forest canopy is difficult to obtain with optical stereo pairs, InSAR and Pol-InSAR technique as the surface is covered with forest and the images hardly include terrain information.DEM beneath forest canopy is an important part of terrain measurement.Both of the forest's structure and the underlying DEM are essential for forest biodiversity and forest ecosystem succession processes.In a broad sense, they belong to forest vertical structure.Monitoring of vertical structures in forest and extracting DEM information beneath forest canopy is required for global carbon cycle characterization and for terrain precise measurement.
Synthetic Aperture Radar (SAR) is a kind of modern high resolution microwave imaging radar.With its long wavelength compared with optical sensors, SAR has the potential to retrieve surface parameters (such as forest map, land cover, target detection and so on) (Toan et al. 1992, Colesanti et al. 2003) due to its easy penetration into vegetation canopy.Polarimetric SAR (PolSAR) facilitates us in investigating and discriminating the scattering mechanisms in forest (Doulgeris, Anfinsen andEltoft 2008, Yamaguchi, Yajima andYamada 2006).SAR Interferometry (InSAR) measures the phase differences between two or more SAR images, which are primarily related to the surface topographic elevation and deformation (Hooper 2008, Wright, Parsons andLu 2004).Polarimetric SAR Interferometry (Pol-InSAR) merges both the abilities of PolSAR and InSAR, in other words, it obtains the scattering mechanisms and phase differences (Garestier et al. 2006, Lopez-Sanchez, Hajnsek andBallester-Berman 2012).SAR tomography (TomoSAR) extends conventional 2-D SAR imaging into a 3-D reconstruction by multi-baseline InSAR / Pol-InSAR images in elevation direction (Reigber andMoreira 2000, Zhu andBamler 2014).Multi-baseline acquisitions are used to reconstruct the reflectivity profile in the direction of elevation, which is perpendicular to the slant and azimuth direction.It leads to a refined analysis of volume structures and terrain beneath forest canopy, and many excellent works have appeared in literatures to explore vertical structures of targets, such as city (Shahzad and Zhu 2015, Wang, Xiang and Wang 2016, Liao et al. 2013, Zhu and Bamler 2014) and forested areas (Li et al. 2016, Bes et al. 2017, Huang, Ferro-Famil and Reigber 2012).TomoSAR, with its several SAR images in different view angles, and sensitivity to semitransparent volumes, has been applied to imaging forest(Stefano Tebaldini 2011) as well as separating scattering mechanisms in layover from urban areas (Wang et al. 2015), localizing targets under foliage (Huang et al. 2012), and imaging subsurface ice structure (Banda, Dall and Tebaldini 2016).
The penetration and different view angles of TomoSAR make extracting underlying DEM possible.In our study, we present an approach for DEM extraction beneath forest canopy by repeat-pass L-band Pol-InSAR data, which is collected in the Krycklan River catchment in northern Sweden in the framework of the European Space Agency (ESA) experiment BioSAR 2008.The paper is structured as follows: part 2, test site and data is described, part 3, methodology for DEM extraction is presented, part 4, results and evaluations are provided, finally, the conclusions and problems are summarized.

Test site
The test site was a forested area in the Vindeln municipality, situated in the province of Vä sterbotten in northern Sweden.The dominating forest type is mixed coniferous forest, the dominating soil is moraine with variations in thickness and the bedrock consists almost entirely of gneiss.The terrain is hilly, with elevation ranges from about 150 to 380 m above the sea level throughout the whole scene.And our test site is considered to be in between a coastal and an inland climate.

Data
The Pol-InSAR data sets have been collected by illuminating the forested areas in the test site at L-band, with look directions South West (SW) and North East (NE), and in this paper only L-band with SW look directions are used.The data set was acquired in the repeat-pass configuration and 6 Pol-InSAR images (Single Look Complex SLC) have been obtained.The campaign has been carried out on October 15, 2008.And the horizontal baseline spacing is approximately 6 m, resulting in a maximum horizontal baseline of approximately 30 m.The spatial resolution is about 2 m in slant range and 1 m in the azimuth direction.The six fully polarimetric SLC images were collected by repeat pass mode, and between the collection no rain or other precipitation.
As a part of the BioSAR 2008 campaign, light detection and ranging (Lidar) measurements were performed by the Swedish Defense Research Agency (FOI).And the DEM extracted by Lidar is applied as the reference data for the validation of retrieved underlying DEM.

Polarimetric covariance matrix
Suppose the geometrics of airborne multi-baseline satisfy uniform or nearly uniform arrangement in cross-range direction, then focused single look complex (SLC) SAR image is obtained by Fourier transform from scenario complex reflectance in cross-range direction, as shown in formula (1).

( , )
( , , ) exp( ) Where, ( , ) Under the natural condition, a simple method is to suppose randomly distributed target changes from one pixel to other, and the data are simplified as realization of the second order moment of the complex process with zero mean.To realize the process, the received signal is supposed as contribution of K different scattering mechanisms, as shown in formula ( 2 Without loss of generality, suppose polarimetric signature scale factor 2 (1) 1, Where, ( , ) k nm  represents the interferometric coherence of w polarimetric channel data collected by the nth and mth track.When the number of polarimetric channels is 3, the covariance matrix polarimetric interferometric covariance matrix, and the polarimetric vector is shown in formula (4). Where, y w y w y w y w  is the projection vector in polarization basis such as Lexicographic or Pauli basis in polarimetric channel 1 2 3 , , w w w .The polarimetric steering vector and steering matrix are shown in formula (5).
A polarimetric steering vector is defined by 5 parameters, including the height z, real and imaginary parts from two complex data.The multi-baseline polarimetric interferometric signals can be represented as: Where, 1 () is the complex amplitude of the ith signal source, i k is the polarimetric target vector of the ith signal source , and the multi-baseline polarimetric interferometric covariance matrix represents as [14].
Where,   is the signal covariance matrix.In practice covariance matrix estimation usually uses effective multi-looked interferometric coherences as shown in formula (8).Sample covariance matrix ˆFP R is able to replace covariance matrix FP R under the assumption that the reflectivity is evenly distributed and sufficient number of looks.
It is important to note that covariance matrix must be normalized, and then all the elements at main diagonal are 1, and FP R is the full baseline interferometric coherence matrix.

TomoSAR
TomoSAR is the extension of conventional 2-D SAR imaging to 3-D.And it is achieved by the formation of an additional synthetic aperture in cross-range (Huang et al. 2012).Usually, this is based on coherent combination of InSAR or Pol-InSAR images acquired from repeat-pass flight tracks of airborne or satellite borne SAR system.TomoSAR directly retrieves the distribution of the backscattered power in the vertical height, leading to section the volume structure, and it is widely applied in forest structure retrieval, building height estimation and so on.
After demodulation, deramping and FFT, the received signal of a scaterer at the position   0 0 0 ,, s r z can be represented by , , exp sin Where, 0 z is the height of the scatterer,   0 0 0 ,, a s r z is the complex reflectivity of the scatterer at the position   0 0 0 ,, s r z , and tomo L represents tomographic aperture and it is the largest baseline length in case of repeat-pass configuration.What should be noted is that the tomographic resolution depends on the width of the sinc-function and it told us a better resolution requires a larger tomographic aperture tomo L .

Ground scattering characteristics extraction
Sum of Kronecker products (SKP) is applied to the sample covariance matrix and only the first 2 terms are retained to identify ground and volume contributions.Then algebraic synthesis (AS) technique is used to retrieve the ground-only contributions by maximizing the corresponding coherence amplitudes in all interferograms under the constraint of physical validity (Tebaldini and Rocca 2012).Usually, ground scattering is supposed to be associated with relative precise location in the vertical height.Then the ground phase in the nmth interferogram is written as () Where, zn k is the wavenumber for the nth track, g z is the ground height and n  is the contribution of propagation disturbances in the nth track.The estimation of 6 phases is resolved by phase linking algorithm and then ground phase is obtained through formula (10).
The estimated ground phases then are processed with standard multi-baseline interferometric techniques to extract phase contributions due to terrain elevation.

Underlying DEM estimation
Phase unwrapping should be performed before extract underlying DEM.As phase unwrapping is not the main concern of the paper, we just used the classic Goldstein branch-cut algorithm to unwrap the phase (Goldstein, Zebker and Werner 2016).Then the unwrapped phase can be applied to estimate the underlying DEM with formula ( 11) and ( 12).
Where,  is the incident angle of the master image, R is the slant range between the sensor and the terrain,  is the intersection angle between the horizontal direction and the baseline, B is the baseline range, H is the altitude of the sensor and  is the wavelength.

Vertical reflectivity profiles
Vertical reflectivity profiles at azimuth direction cut in HH/HV/VV polarimetric channel are extracted by TomoSAR technique.And the profile position and vertical reflectivity profiles are shown in Figure 1.We could see from figure 1(b) that the vertical reflectivity varies in the height and the distribution of reflectivity is in line with the in-situ target.The profiles in dense forest range at higher location and the height of profiles in secondary forest are near the ground.It is noted that the vertical backscattering power in each polarimetric channel seems to rise from near to far range.This phenomenon still exists, although MUSIC technique has been applied in vertical backscattered power reconstruction.This may be caused by the large variations of the normal baseline span along the imaged swath, and then results in the vertical resolution to vary as well.

Ground and volume structures
Using SKP and AS technique, we obtained ground and volume scattering characteristics, also called structure.As scattering characteristics separation is based on algebraic operation and no physical scattering mechanism is considered.There is mutation between the extracted ground and volume scattering characteristics, as shown in figure 2 (a-d).However, it is more important is to identify ground scatterers in this part.And the figures indicate that the reflectivity ranges of extracted ground structure are near the ground, and it is consistent with the theory of the distribution of the classical surface scatterers.

CONCLUSIONS AND DISCUSSIONS
We has presented evaluation of L-band TomoSAR for underlying DEM estimation in the boreal forest within the Kryclan test site, investigated in the frame of the ESA BioSAR 2008 campaign.It has been observed that the extracted ground structure could reflect the trend of terrain, which is consistent with scattering mechanism separation using AS algorithm in polarimetric TomoSAR (Tebaldini and Rocca 2012).The general trend of estimated DEM beneath forest canopy is the same with that of Lidar DEM, and the mean difference is about 3.39 m for our test site.Some of the estimated DEM is lower than Lidar DEM, which may relate with that the penetration of L-band is better than that of laser in forested area.The work indicates that L-band polarimetric TomoSAR is feasible for underlying DEM extraction in forested areas.
TomoSAR represents a powerful tool for investigating the distributions of reflectivity and SKP+AS algorithm provide a useful way to separate scattering characteristics in vertical direction.Both of these make underlying DEM estimation with L-band repeat-pass Pol-InSAR data possible.However, there are some problems to be studied further.First, the physical target decomposition will be considered in scattering mechanism separation to enhance separation degree.Second, the approach of phase extracted from the ground structure needed to be optimized further and P-band data will be used to decrease the effect of forest canopy.
n y r x represents the pixel in the nth SLC image, ( , , ) s r x v is the mean complex reflectance of 2D SAR image at   , rx , n b is the spatial baseline of the nth image, and v is the resolution in cross-range.
To sum up, the covariance matrix for single polarimetric TomoSAR is shown in formula (3). 1 (a) Pauli-bases image and position of profiles (Red line represents profile position at slant range, and blue line represents profile position at azimuth direction) (b) Vertical profiles in HH/HV/VV polarimetric channel at slant range Figure 1.Positions and vertical reflectivity profiles (a) Ground structure at azimuth direction (b) Volume structure at azimuth direction (c) Ground structure at slant range (d)Volume structure at slant range Figure 2. Ground and volume structures at profile positions in figure 1 (a) 4.3 Underlying DEM The underlying DEM is estimated from the separated ground structure and the permanent scatterers extracted from intensity image.The extracted underlying DEM and Lidar DEM are shown in figure 3. From the figure, we can see that the range of the underlying DEM is 200 and 280 m, while the Lidar DEM is from 190 290 m.The general trend of elevation of the two DEM is the same, although there are some different details in the image.Figure 4 (a) is the distribution of differences and (b) is the according Pauli-base image.From the figures, we could see that maximum difference locates in the complex terrain (large slope areas and dense forest).The mean difference between the estimated underlying DEM and Lidar DEM is 3.39 m, and maximum difference is about 30 m. Differences between Lidar DEM and estimated DEM(a) and Test site Pauli-base image (b)