TERRAIN CORRECTION OF POLARIMETRIC SAR DATA AND ITS APPLICATION IN MAPPING MOUNTAIN GLACIER FACIES WITH RADARSAT-2 SAR

Synthetic Aperture Radar (SAR) data can be difficult to interpret over glacierized terrain, especially in areas of high relief where most of the smaller glaciers of the world are found. The two primary objectives of this paper are to correct the topographically induce distortions inherent in SAR image of rugged terrain and to map glacier facies in the corrected polarimetric (PolSAR) image. Methods for terrain correction of PolSAR developed in this paper, which base on the coherency matrix data, utilize polarimetric signature to refine the location accuracy and perform the pixel size normalization on each element of the coherency matrix using the backward integration method. Moreover, azimuth-slope correction is immediately following after radiometric correction. A supervised classification is performed on the orthorectified image to map the spatial distribution of snow and ice. The visual boundary between wet snow and bare ice surface is delineated on the orthorectified image too. These results show that glacier facies and the wet snow line can be mapped reasonably well using Radarsat-2 PolSAR in a mountainous environment. * Corresponding author.


INTRODUCTION
As it not affected by cloud covers and has ability to penetrate beneath to ice/snow surface to observe the subsurface structure, SAR has been successfully applied to monitor the Mountain Glaciers due to its high spatial resolution and wide coverage (Henderson and Lewis, 1998).However, spatial and temporal variations in snow pack make SAR backscattering signal for glacier facies poorly understand, including the relative important of scattering from surface and volume layers (Pellika and Rees, 2009).Recently, fully PolSAR systems such as Radarsat-2(C-band) and PALSAR (L-band) have been launched.These new instruments promote the detection and monitoring capabilities of glacier facies as the polarimetry contain information on the structural characteristics (shape and orientation) of the scatterer and enable us to better discriminate between volume and surface scattering and furthermore to properly indentify the different glacier facies.
However, when trying to monitor the snow characteristics and map glacier facies from PolSAR data in a mountainous environment like Mountain glaciers in western China, radiometric corrections must first be applied to correct for the distortions induced by the slant projection of SAR systems and by the highly variable terrain.The terrain effect can be split into two components: the change of radar cross section per unit image area and the polarization state since the azimuth slopes induce polarization orientation changes (Lee et al., 2000).
The objective of this paper is twofold: we want to extend the theory of SAR terrain correction to polarimetric case, where correction is performed explicitly based on matrix format as multi-looked PolSAR data are usually in formats like covariance, coherency, or Stokes matrix, and therefore utilize the all polarimetric signature available.We next want to apply the ortho-rectified PolSAR to map the glacier facies and investigate the potential of PolSAR in mapping glacier facies.First, the radiometric slope correction is performed by using the rigorous SAR image simulation and mathematical modelling of SAR geometry.Specially, the all polarimetric signature is utilized to refine the location accuracy and then these projection factors are adopted to normalize the all elements of coherency matrix in map space using the backward integration method, which taking into account the slope and aspect of the DEM and preserving the energy contained in the image data throughout the geo-coding process.Secondly, each element of the PolSAR matrix is compensated according to its polarization state by rotating the data along the line-of-sight by the negative of the induced orientation angles.Finally, a supervised classification is performed on the orthorectified PolSAR to map the spatial distribution of snow and glacial ice and evaluate the detection and monitoring capabilities of glacier facies with PolSAR.
In the following sections, we first introduce the test site and experimental datasets, and then describe the algorithms and processing steps involved in our method on terrain correction and mapping glacier facies, and finally, present some technical remarks and conclusions.

TEST SITE AND EXPERIMENTAL DATA
The test site is the Dongkemadi Glacier Basin located in the in the middle of the Thang-la Mountain Range of Tibetan Plateau of China, which covers an area of 50.96km 2 ranged in elevation from 5000m to 6100m above sea level (a.s.l.).The annual average temperature is -6 ℃ and only from June to September is above 0 ℃ in this drainage.In winter, under the westerly circulation control, the dry sunny and windy climate can last for 8 months (from October to May of the following year); In autumn, the warm and moist air blow from Southwest Indian and the precipitation concentrate in the range June to September (Fujita et al., 2000;Shangguan et al., 2008).The valley formed by the ancient glacier is flat and open.The surface of glacier is flat, clean and no debris.Meltwater channels in the bare ice facies appeare only in late summer when the strongest melting happens.To our known, the retreat state keeps on since the summer of 1994 and it had retreated 4.56m by 2001. Figure 1 show the satellite images of Dongkemadi Basin.The Dongkemadi Glacier locates in the red rectangle of figure 1 A digital elevation model (DEM) of Shuttle Radar Topographic Mission (SRTM) with 90m pixel spacing covered Dongkemadi Basin was downloaded from the website of USGS U.S.A.In order to maintain a certain spatial resolution, the DEM data was oversampled in range and azimuth using the oversampling factors 2. The DEM data mainly is used for orthorectification and generate contour line for analysis.

TERRAIN CORRECTION
Because of the side looking geometry of a SAR system, every target located on the terrain being observed by the radar will be mapped onto the slant range domain (Schreier, 1993).Radar imagery is likely to be affected by geometric distortions and the most common distortions are these of foreshortening, layover and shadow.The topographic effects caused by rugged terrain have two components, namely a geometric and a radiometric one.In geometry, the topographic effects result the relationship between radar and map geometries is not homographic.
Foreshortening and layover induce a many-to-one mapping, namely, many surface target points may be mapped to one pixel in the slant range image and ambiguity in the reverse direction also occurs in proximal slope.In addition, the scattering area for each pixel may be a complex and some noncontiguous area of ground and thus accurate radiometric correction is essential (Luckman, 1998).Using a DEM, the imaging geometry is reconstructed and then is used to perform geometric and radiometric correction of terrain induced distortions.

Imaging geometrical model
The objective of the accurate reconstruction of SAR imaging geometry is to find for each image pixel, the corresponding position on the earth.In an earth centred Cartesian coordinate system, for given a target position P , we need to determine its two-dimensional image coordinates: ( , )  ij .In slant range geometry, the image row number i is known as the azimuth (along-track) coordinate, and the image column number j as the slant range (cross-track) coordinate.With the known orbit data vector and image-line timing, the geocoding lookup tables (LUTs) can be obtained by solving the range and Doppler equations (1) for each DEM cell (Curlander and McDonough, 1991).
where R is the slant range from the SAR to the target; P is the known position of target in 3D Cartesian coordinate system; t is the time when the target is imaged; S is the SAR position at time t; d f is the Doppler frequency shift of the target;  is the radar wavelength; P V is the velocity vector of the target due to the earth's rotation; S V is the velocity of the SAR; and the ()  denotes the inner product of vectors.The solution of this nonlinear system is found through iteration along the orbit for each DEM-pixel with the Newton-Raphson method.The results of the iteration (imaging time and range coordinates) are linearly transformed into coordinate system of the slant range image and thus obtain the initial geocoding LUTs.In addition, the reconstruction of imaging conditions is done to derive various geometrical parameters for each ground point (e.g. the look angle, the local incidence angle and the projection angle etc.) and generate the layover-shadows map (Ulander, 1996).Figure 2  The initial LUTs are not completely accurate because of limitations in the knowledge of the SAR imaging geometry (the orbit positions are only known with certain accuracy).As a result geocoding with the initial lookup table leads to geocoding errors of the order of a few pixels, depending mainly on the accuracy of the orbit data.In order to refine the lookup table, the simulated SAR image is applied to indentify Ground Control Points (GCPs) by image matching, which is performed between the simulated image and the real SAR image.Based on these obtained GCPs, an affine transformation model is constructed to characterize the mapping relationship between the simulated and real SAR image coordinate space.It is different from traditional methods that we make a polarimetric image fusion and generate a composite image, which may provide useful help for real-simulated image matching (e.g.preserve image contrast and enhance texture while remove the speckle).(Toutin et al., 2010) have evaluated four different kind of polarimetric image fusion algorithms by a variety of image quality measures (equivalent number of look (ENL), texture, mutual information and auto-and cross-correlation coefficients), and concluded that three fusion algorithms (SPAN, PC analysis and restoration) were considered to be the most appropriate for image matching.However, in order to utilize the all elements of T 3 and thus all polarimetric information, we selected the mean eignvalues  calculated by Cloud-Pottier's decomposition method (Cloude and Pottier, 1996) as the fusion image.For quantitative estimation of the geocoding accuracy, manual selection of independent check points (ICPs) is interactively performed with the aid of a graphic tool in the ENVI image environment.Figure 3 shows the 30 ICPs vectors exaggerated by a factor 100, which are the differences between the known and the calculated image coordinates for each check point.A least-squares adjustment with affine geometric mode resulted in RMS of 0.667 pixels in azimuth direction and 0.59 pixels in range direction, which is better than the pixel spacing of SAR image.

Radiometric Correction
Previous researches have shown that for SAR imagery from noflat terrain, the radiometric correction must take account the pixel scattering area, which is dependent on the local topography as well as the look angle of the radar.After an accurate reconstruction of the local SAR geometry, these parameters can be used for a precise radiometric correction of the SAR image.Due to the RADARSAT-2 have a better spatial resolution (about 8m) than the SRTM DEM (90m), for rugged glacier area the scattering area corresponding to each image pixel can't be estimated accurately in slant range space and select an output spatial resolution much higher than the input spatial resolution strongly reduces the geocoding efficiency without much information gain.Therefore, we first apply the backward integration method to calculate the radar brightness for each DEM pixel, then perform the radiometric normalization in map geometry and finally compensate the PolSAR data for terrain azimuth slope variation.

Radiometric Normalization:
To calculate the corresponding radar brightness for each DEM pixel, the traditional backward resample method assign the image pixels to the DEM coordinates by using bilinear or cubic convolution kernels.In the case of SAR image with a better spatial resolution than the DEM, due to the sampling is inadequate for the interpolator, the signal noise level is kept unnecessarily high in resampling SAR images from a higher resolution to a lower resolution.To reduce speckle effect when reduce the spatial resolution to a similar value as required for the output image (e.g.DEM resolution), the integrate method is applied in this paper.We first assign all image pixels to its corresponding DEM pixel according to the closest Zero-Doppler position and then averaged.For single-polarization radars, the average radar brightness can be compute: where all N image pixels share the same Zero-Doppler coordinate on the ground and each image pixel is assigned once to a DEM pixel only (Loew and Mauser, 2007).For a PolSAR, the integration can be made in an analogous way using T 3 .This method guarantees that the integrated backscatter intensity, measured for each pulse, is preserved throughout the geocoding process which is essential for the generation of geocoded products, which are comparable to the original SAR image and is the basic requirement for a successful radiometric normalization of the backscattering coefficient.In order to quantify the influence of the integrate method, the frequency distributions of the T 11 elements of T 3 in glacier region were calculated with two different method and the result are shown in Figure 4. We can find that there are significant difference frequency distribution and the average value decrease from 1.31 to 1.27 (the standard deviation from 1.88 to 1.36) due to integrate.This corresponds to the large part of the investigate area with height slops, in where the backscatter intensity measured for each plus have a very large difference.The average process reduces the impact of this dissimilarity and radar speckle.The normalized radar backscattering coefficient is related the radar brightness (Ulander, 1996) as 00 σ β cosψ  (3) where the projection angle is intersection angle between the surface normal of the slant range plane and the ground surface and cosine projection angle , namely, the projection factor can be derived from SAR look vector and local terrain slops and aspects as cos ψ sin θ cos u cos θ sin u sin v  (4) where is the local incidence angle, u and v are the terrain slope and aspect with respect to the sensor azimuth direction separately.The projection angles are shown in Figure 5(a) and its corresponding projection factors shown in Figure 5(b).We can see that the backscattering intensity of face the radar is reduced when using the projection factor to normalize the coherency matrix.Moreover, for the projection angle lager than 90 degree, the projection factor is zero, corresponding to the deepest blue areas in Figure 5(b).These areas are affects by layover and have to be masked during geocoding process.

Azimuth slope correction:
Lee et al have point out that the azimuth slopes affect the relative magnitude and phase of all terms of the polarimetric coherency matrix.Thus, for terrain classification and other PolSAR application, the azimuth slope effect should be corrected to ensure accurate estimation of geophysical parameters (Lee et al., 2000).The polarization angle change due to the azimuth slope effect can express as tan ω tan tan γ cos sin where is the azimuth slope, is the slope in the ground range direction, and is the radar look angle.The polarization orientation shift is not affected not only by azimuth slop but also by the range slop.After derived the orientation angles (OAs), the compensation can be done by 33 new T  T UT U (6)   where the rotation matrix U is given by 1 0 0

U
Lee recommended that azimuth slope compensation be performed on L-band and P-band data not on C-band and higher frequency data.They consider that the orientation angles associated with small scatterers overwhelm that from the ground slop in C-band.However, the most regions of our test site belong to glacier without vegetation cover.By extracting a profile of the polarization orientation angels along the central flow line of Big Dongkemadi Glacier (from its accumulation area down to its ablation area, see the yellow line on B in figure 10), we investigate the difference of the orientation angles derived by DEM directly and circular polarization method.The results are shown in figure 6.For circular polarization-derived OAs, the noise effect is quite visible although the trend is consistent with the DEM-derived OAS.Due to the wavelength of C-band radar is close to the scale of surface roughness, the scattered signal is more sensitive to small-scale roughness and cannot remain correlated over the correlation length of the large-scale roughness.So, we use the DEM-derived OAs to compensate the terrain azimuth slope variation in this study.

Results in terrain correction
The geocoding accuracy was evaluated by visual checking the overlap between DEM and the SAR image in figure 7.Although no GCPs were used to measure its accuracy exactly in a quantitative way, we are able to find that the lines of ridge and valley are agreed very well and these traits of radar geometry are highly correlated with the topography, e.g.layover.

MAPPING GLACIER FACIES
As one can see in figure 1, by presenting the Pauli-based decompositions in RGB, various land cover objects are emphasized.The stone crushed by glacial are clearly revealed in green due to its extremely rough surface and strong depolarization caused by multiple scattering.Although sparse and short (about 3-5cm) grasses grow on the flat valley, surface scattering mechanism plays the dominate role, which can be judged from the blue tone of these region.In the wetsnow facies, the radar signal is absorbed by the thin layer of liquid water, resulting in a dark (low) return, while in the bare ice facies, rough surface provides a stronger backscattering.Therefore, the boundary between the wet-snow facies and the bare ice facies is detectable on radar image.In the following section, we will investigate the capability of Polarimetric Radarsat-2 to map the glacier facies by using the image after terrain correction.

Supervised classification
The basic idea of SVM is to map data into a high dimensional space and to find the hyperplane the different classes with the maximal margin between them.The attractiveness of SVMs is their ability to minimize the so-called structural risk, or classification errors, when solving the classification.Assume there are two linearly separable classes (Tso and Mather, 2009).
The training data set is represented by pairs   where training data are mapped to a higher dimensional space by the function  , and C is a penalty parameter on the training error.Practically, we need only , the kernel function, to train the SVM.In this paper, the kernel function was applied is the radial basis function (RBF): Polarimetric SAR is becoming more widely used in remote sensing classification, especially; many kinds of polarimetric parameters were applied to classification (Cloude and Pottier, 1997;Lee et al., 2004).There is an obvious link between the number of input feature and the computational burden, yet classification accuracy may not show any significant improvement if one of the input features dose not significantly contribution to the classification problem.After an in depth research on the suitability and the usability of the parameters proposed in the literature, we selected all coefficients of Paulibase, entropy and anisotropy as the PolSAR features to classify the research region.
Figure 9. results of classification by using SVM.

Find the wet snow line
Many studies on smaller glaciers find a distinct boundary between wet snow and bare ice on the SAR image at the end of a mass-balance year, the boundary is either transient snowline or the firn line, depending on the weather conditions during the ablation season (Jaenicke et al., 2006;Konig et al., 2001).According to our field investigation, only the wet snow facies and the bare ice facies can be observed on Dongkemadi Glacier and no firn exposed.Thus the position of wet snow line contains information about the mass balance as the equilibrium line equals the wet snow line approximately at the end of a balance year.
First, these positions of wet snow lines were investigated by overlaying the high contours in figure 10.We can find that the snow line altitude is range from 5550m(blue line in figure 10) to 5600m(red line in figure 10) (a.l.s.). Figure 10 also shows that the wet snow line is at a higher altitude on the South (A and B) than on the North (D, E and F).This is probably due to higher direct in indirect sun irradiance because of its southfacing slopes.Secondly, in order to find the wet snow line more accuracy, the Optimization of the Polarimetric Contrast Enhancement (OPEC), which enhance the contrast between two types of scatterers by finding the optimum polarizations, was performed on bare ice (target) and wet snow (clutter) areas of interest (AOI).After image enhance, the snow lines were delineated manually on six outlets of glacier.Figure 11 illustrates the snow line altitude on various outlet of Dongkemadi Glacier.

DISCUSSION AND CONCLUSION
This study shows a rigorous method for geometric and radiometric correction of Polarimetric SAR data.There are three improvements in terrain correction: (1) refine the location accuracy by polarimetric image fusion; (2) extend the integration method of the radar brightness to polarimetric case for radiometric terrain correction; (3) investigate the azimuth slop compensation for Radarsat-2 on glacier region.After terrain correction, we applied the corrected PolSAR data to map the glacier facie with supervised classification SVM and illustrate the potential of RADARSAT-2 in monitoring the Mountain glacier in Tibetan of china.

ACKNOWLOGEMENTS
This work was supported by the Chinese Ministry of Science and Technology (Grant Numbers 2009CB723901, 2010CB951403 and 2006FY110200) images of the Donglemadi Basin.(a)Paulibaesed colored PolSAR image; (b) the optical image (from Google Earth) The satellite imagery used for this study is RADARSAT-2 Cband PolSAR data, which acquired on 30-Aug-2009 at 11:48 UTC time.The pixel size of the data is 4.7m in range direction and 5.5m in azimuth direction.The imagery has undergone relative calibration at Canada Centre for Remote Sensing (CCRS) prior to distribution.The first processing step for SAR images is converting digital values to radar brightness β 0 and then saving in the 3*3 coherency matrix T 3 format.

Figure 2 .
Figure 2. (a) DEM of test site; (b) the local incidence angles.
(a) shows the DEM of covered area and figure 2 (b) gives the local incidence angle.

Figure 4 .
Figure 4. Histogram of element T 11 of coherency matrix calculated by (A) backward resample and (B) integrate method.

Figure 5
Figure 5 the radiometric normalization parameters of the covered area.(a) for the projection angles and (b) for the projection factors.
Figure 6.OAS derived by circular polarization and DEM method along the central flow line of Big Dongkemadi Glacier.

Figure 7 .Figure 8 .
Figure 7. 3D DEM draped PolSAR image after correction To visualize the effect of data terrain correction, line profiles form the area B in figure 10 are plot in figure 8 for the diagonal elements of coherency matrix, before and after correction.The original radar brightness values are shown in red broken-line, the radiometric normalized values are shown in green line, the compensated with DEM values are shown in blue dotted-broken line and the after compensated with circular polarization values are shown in black dotted line.The power of |HH -VV| in figure 8(a) is consistently increased after the OA compensation, while the power of HV shown in figure 8(b) is correspondingly reduced and the |HH + VV| is roll invariant in figure 8(c).This confirm the observations made by the previously cited authors

Figure 10 .
Figure 10.PolSAR image overlapped high contours.(A-Fdenote the location of six outlets of glacier analyzed)