DETERMINATION OF GLACIER SURFACE AREA USING SPACEBORNE SAR IMAGERY

Glaciers are very important climate indicators. Although visible remote sensing techniques can be used to extract glacier variations effectively and accurately, the necessary data are depending on good weather conditions. In this paper, a method for determination of glacier surface area using multi-temporal and multi-angle high resolution TerraSAR-X data sets is presented. We reduce the “data holes” in the SAR scenes affected by radar shadowing and specular backscattering of smooth ice surfaces by combining the two complementary different imaging geometries (from ascending and descending satellite tracks). Then, a set of suitable features is derived from the intensity image, the texture information generated based on the gray level co-occurrence matrix (GLCM), glacier velocity estimated by speckle tracking, and the interferometric coherence map. Furthermore, the features are selected by 10-foldcross-validation based on the feature relevance importance on classification accuracy using a Random Forests (RF) classifier. With these most relevant features, the glacier surface is discriminated from the background by RF classification in order to calculate the corresponding surface area.


INTRODUCTION
The glaciers in mountain areas sensitively interact with climate fluctuations and therefore they become very important climate indicators.Satellite remote sensing has a great capability in monitoring the glaciers located in cold high altitude regions and inaccessible areas.Visible sensors have been extensively used to map ice covered areas (Paul et al., 2002;Paul et al., 2004;Kargel et al., 2005;Hendriks and Pellikka, 2007), but the limitations are their dependence on weather conditions, which cannot provide timely information of the target area, especially some mountainous glaciers with year round cloudy conditions.However, SAR, due to its all day, all weather imaging capabilities, can reliably collect data with a pre-defined temporal interval over long periods of time with a spatial resolution compatible with the application of glacier surface area determination.Within opposition to that, using optical remote sensing imagery it might happen that cloud-coverage makes one (or more) acquisitions in the interval impossible, thus leading to a time series with "gaps".Much work has been done to map and study glaciers and ice sheets using SAR techniques: most of the work is based on interferometric information or polarimetric data (Shi and Dozier, 1993;Dall et al., 2004;Zhou et al., 2010).In this paper, we propose a method for glacier area determination using multitemporal and multi-angle high resolution TerraSAR-X data sets aiming to reduce the 'data holes' caused by radar shadowing and specular backscattering of smooth ground surface.The surface area of a glacier is an important parameter not only for the two-dimensional mapping of changes of the glacier extent; it furthermore can be used to get rough estimates of the corresponding glacier volume based on geophysical models (Bahr et al., 1997).

GLACIER SURFACE ESTIMATION
The work proposed in this paper aims at the estimation of glaciers surface areas.For this, first a supervised classification is carried out in order to distinguish the glacier from the background.Afterwards, the classification result is filtered and post-processed using low-level morphological operators.Based on this procedure, the glacier outline is delineated and the surface area calculated.The glacier velocity, an important control parameter determining the mass balance, is a crucial feature that discriminates the glacier from the surrounding, non-moving, environment.The glacier motion maps utilized for the work described in this paper are obtained by speckle tracking based on cross-correlation of the intensity images.However, many 'data holes' appear in these motion if they are generated from SAR pairs imaged from just ascending or just descending satellite tracks: On the one hand, mountainous terrain is affected by radar shadowing, on the other hand smooth ice surfaces cause specular backscattering such that no signal is returned.Thus, based on multi-temporal and multi-angle SAR images (at least 4 scenes with 2 from ascending track and 2 from descending track), the motion maps derived from both imaging geometries are fused using a probabilistic graphical model for improving the estimation coverage and accuracy.The possibility to include context knowledge into the probabilistic graphical model furthermore helps to make the results more reasonable concerning ice physics.The potential for the improvement of ice motion estimation has already been shown in (Fang et al., 2013).Figure 1 shows the final velocity estimation result from a part of the surging area of the Taku Glacier.It is supposed that the velocity information should be a very good indicator for discriminating the glacier from the nonmoving scene content (e.g., mountain area).

Coherence Map:
The interferometric coherence map is generated by using two TerraSAR-X repeat cycle scenes acquired in same imaging geometries (either both ascending or both descending).As the moving glacier shows strong temporal decorrelation leading to low coherence values, the coherence map can be an important indicator for glacier classification.

Texture Information:
Texture is one of the important characteristics used in image classification.The second-order gray-level co-occurrence matrix (GLCM) proposed by Haralick (1979) is a well proved effective approach to analyzing image texture features.The GLCM approach to texture analysis is based on the conjecture that the texture information in an image is contained in the overall or average spatial relationship between the gray tones of the image.In this work, the GLCM of each pixel is generated by 4 directions comparing its 14 neighbourhoods which surround this centre pixel.Afterwards, 4 features are derived from the GLCM: contrast, correlation, energy and homogeneity, which are calculated by Then, also the mean and standard deviation values of the 4 calculations from 4 different directions of these 4 features are calculated, such that finally 8 texture features are defined to be used in the experiments.Much of the glacier surface area is featured by characteristic textures (e.g., crevasses or drainage patterns), which makes the texture a useful kind of information for the classification.

Others:
The other basic features used in this investigation can be directly retrieved from the image data: the intensity value and the other 3 features derived from the intensity value employing a 7*7 window size (mean, median, standard deviation).

Random Forest:
The Random Forest (RF) is an ensemble learning method for classification (and regression) that is based on a multitude of decision trees at training time and outputting the classification result that is the mode of the classes provided by individual trees.Each tree in the forest classifies the sample based on an independently sampled random variable set which are used to best split the node under a certain criteria (e.g., information gain, Gini impurity).After every tree grows until minimum node size is reached, the forest chooses the classification having the most votes (Breiman, 2001).In parallel, the feature relevance (variable importance), i.e. the contribution of the respective feature to the classification accuracy is generated at the training time, which can be a useful measure for feature selection and reduction of the feature space dimension.

Feature Selection Method:
This work utilizes a feature selection method based on feature relevance generated by RF and the sequential backward selection method (Wei et al., 2012).The general idea is that we rank the features firstly according to the feature importance score while the overall classification accuracy is calculated for this variable set, and then delete the feature that scored lowest.This operation is repeated several times until we find the variable set with the best classification performance.In order to make the test more reliable, in every round for feature selection, 10-fold-cross-validation is used.The samples are randomly divided into 10 parts, with nine parts for training the tree classifier and the rest for test.During the 10 loops, the relevance importance in the loop with highest classification accuracy for the test samples is used for ranking features.Additionally, the average classification accuracy of this round is the mean value of the 10 loops.

Glacier Classification
Combining the selected features based on the method described above, a classification of the glacier surface can be carried out by the RF classifier.As the RF classification constituted by two parts (i.e., training and testing parts), we use two separate SAR scenes, with the same area coverage, for training trees and final classification, which on the other hand can also guarantee the experiment reliability.

Glacier Outline Delineation
Aiming to generate the reliable delineation of the glacier used for an accurate surface area calculation, the classification result is firstly median filtered, followed by the morphological International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-1/W1, ISPRS Hannover Workshop 2013, 21 -24 May 2013, Hannover, Germany operations closing and opening.Finally, the outline is detected by tracing the exterior boundary of the glacier area in the binary classification map.

Glacier Surface Area Calculation
According to the glacier classification result, the glacier surface area surface area A can be roughly estimated easily by multiplying the number of correctly classified glacier pixels N with the pixel size: * * A N a b = (5) where a and b are the pixel size in azimuth and range direction, respectively.A more precise estimation of the glacier surface area can be achieved, based on the glacier outline detection result, using the Gauss method for area calculation: ( ) where x and y are the 2D positions of the points describing the border of the glacier.

STUDY AREA AND DATASETS
The Juneau Icefield is a low-latitude glacier system of small scale located in southeast Alaska, and much research work has been carried out to study the glaciers located in this Icefield (McGee et al., 2007;Pelto et al., 2008).In this paper, the scenes of Juneau Icefield imaged by high resolution TerraSAR-X are applied for experiments.All the experimental works in this paper are based on 4 TerraSAR-X repeat cycle scenes acquired in stripmap mode from two different imaging geometries (ascending and descending).2 scenes from ascending orbit are acquired with 11days interval, and the other 2 imaged from the descending orbit are complementarily used for reducing the "data holes" in the SAR scenes (see Table 1).The TerraSAR-X Single Look Slant Range Complex (SSC) data used are firstly orthorectified and geocoded using precise available DEM data.The enhanced data have a ground resolution of 3m with pixel spacing of 1.5m in both azimuth and range direction.The processed SAR scenes used in this research are plotted on the ground truth area shown in Figure 2. Finally, one part of the Taku Glacier surging area marked by a yellow square box is cut aiming to test the method for delineating glacier surface addressed in the paper.For the classification using the RF method, we use one of the ascending SAR images as the training data for tree model construction.Then, the trained tree model will be applied to the other ascending SAR image for classification.
The test scene of the Taku Glacier (marked by the yellow box in Fig. 2) used for classification is shown on the upper part of Fig. 3.The manually generated ground truth mask image is shown on the lower part.All in all, 4941038 pixels belong to the glacier, while 6455851 pixels belong to non-glacier areas.

EXPERIMENT RESULTS
4.1 Generated Feature Set 4.1.1Velocity Image: Using the method described in section 2.1.1,the velocity image in the test area is estimated with a smoothness constraint, reasonably considering the glacier physics (Fig. 4).Furthermore, the motion value differences between glacier and its surrounding neighbourhoods are helpful for classification.

Coherence Map:
The coherence map generated from the ascending SAR pair from the datasets is show in Fig. 5.It was found that the glacier area shows relatively low coherence, appearing as black in the coherence map, whereas the surrounding mountainous area shows higher coherence.However, the water area also shows low coherence, making it a less useful indicator for the classification of glaciers.

Texture Information:
Figure 6 shows the mean values of the four features calculated from GLCM by ( 1)-( 4) as described in section 2.1.3.The water area is recognised as dark area in contrast, while bright in energy and homogeneity, which is due to the fact that water areas don't show signal information in SAR images because of the specular backscattering.This is helpful for discriminating the glacier from the water area.Additionally, the standard deviations of these four calculations from GLCM are shown in Fig. 7 as a complementary measurement.

Feature Selection Result
Based on the method described above, the original 14 features are processed, and finally a set with 11 features achieving 91.05% percent of classification accuracy is selected as the best choice.

Glacier Classification Result
Furthermore, based on the 11 selected features and the RF technique, a classification of the glacier surface is carried out with an overall accuracy of 93.72 percent achieved (Table .2 and Fig. 10).

Glacier Surface Outline
Use the method mentioned in section 2.4, the outline of the glacier surface is generated and shown in figure 11.

Glacier Surface Area
Employing (5), we roughly obtain 10.31km 2 as the first surface area estimation result.Then, according to ( 6) and the surface outline of section 4.4, the more accurate glacier area result is calculated as 11.01km 2 .This result is much different from firstly calculation, which is probably because of the subpixel accuracy achieved by (6).
In order to compare with the ground truth, we plot the four corners of the test scene on optical image via Google earth and then delineate the borders of glacier manually, and the area estimation result is 12.12km 2 (Fig. 12), which has 9.16% relative difference between the result calculated by our method.As the ice balance has been nearly stable for a long period (Pelto et al., 2008), the ground truth data acquired in summer 2010, one year later than our experiment data, is comparable.Then, the relative difference above indicates that the addressed method in this paper for glacier surface determination is effective. International

DISCUSSION AND CONCLUSION
This paper proposes a framework to fuse multi-temporal TerraSAR-X stripmap data from both ascending and descending orbits for the determination of the surface area of glaciers.In order to achieve this goal features such as coherence map, texture information, and pre-estimated velocity maps are used for a classification of the glacier surface.Low-level image processing techniques are then applied in order to delineate the glacier outlines, from which the glacier surface can be estimated.It has been shown that these features allow for an appropriate discrimination between glacier and non-glacier areas.The estimated surface area does effectively fit the ground truth data.However, there are still many improvements left for further research.For instance, the classification of the glacier still fails in some places (cf.red box in Fig. 11), which is because of the lack of velocity and other features in these areas affected by radar shadowing or specular backscattering in both satellite imaging geometries.This might be figured out by detecting glacier edges which are indicated by double bounce backscattering using an appropriate edge detection method and fusion with the classification result.
In future work, a precise knowledge of the glacier surface area is expected to help the mapping of two-dimensional glacier changes as well as serve as an important indicator for glacier volume estimation within geophysical models.

Fig. 1 .
Fig.1.Improved glacier motion estimation result from a part of the surging area of the Taku Glacier.Red dash line: the manually extracted boundaries of the glacier.Although the velocity map already gives a good approximation of the whole glacier surface, there are still some gaps left (e.g. in the upper right part of the image).

Fig. 2 .
Fig.2.Ground coverage of geo-coded SAR data plotted on an optical image of the scene.Red arrow: azimuth direction of ascending orbit; Green arrow: azimuth direction of descending orbit.In both tracks, TerraSAR-X is right looking.Yellow box: area used for experiment.(Underlying optical image ©2013 Google Earth)

International
Fig.4.Fused estimated glacier velocity map plotted on the gray scale SAR image.

Fig. 5 .
Fig.5.Coherence map of one part of the outlet area of the Taku Glacier generated from two TerraSAR-X repeat cycle scenes with 11 days interval.

Fig. 6 .
Fig.6.Upper left: mean value of contrast.Upper right: mean value of correlation.Lower left: mean value of energy.Lower right: mean value of homogeneity.

Fig. 7 .
Fig.7.Upper left: standard deviation of contrast.Upper right: standard deviation of correlation.Lower left: standard deviation of energy.Lower right: standard deviation of homogeneity Fig.8.Relationship between classification accuracy and feature number This eventually leads to the deletion of the 3 least-relevant features (i.e., the intensity value, standard deviation of correlation and energy of GLCM).The importance to the classification accuracy of the remaining 11 features is shown in the Fig. 9.The velocity information is the most important variable, with a rather high score, for classification in comparison to the others, confirming the hypothesis of section 2.1.1.The other features are almost comparably relevant scored far behind that of velocity.When only use the velocity as the feature for classification, the accuracy is 75.73%.

Fig. 10 .
Fig. 10.The glacier surface classification result (red colour) plotted on the geo-coded gray scale SAR scene.

Fig. 11 .
Fig. 11.Glacier outline (yellow line) determination result plotted on the gray scale SAR image.Red box: glacier area failed to be classified.
Fig. 12. Ground truth data plotted on an optical image of the scene.Red polygon: the exact glacier area derived based on the optical image (acquired in summer 2010) manually.(Underlying optical image ©2013 Google Earth)

Table 2 .
Random Forest classification of glacier and non-glacier area