Modelling glacier elevation change from DEM time series

(1)

Master Thesis, Department of Geosciences

Modelling glacier elevation

change from DEM time series

Di Wang

(2)

(3)

Modelling glacier elevation change from DEM time series

Di Wang

Master Thesis in Geosciences

Discipline: Physical Geography, Hydrology and Geomatics Department of Geosciences

Faculty of Mathematics and Natural Sciences

University of Oslo

01.12.2014

(4)

© "Di Wang", 2014

This work is published digitally through DUO – Digitale Utgivelser ved UiO http://www.duo.uio.no

It is also catalogued in BIBSYS (http://www.bibsys.no/english)

(5)

Abstract

Assessment of glacier surface elevation evolution is crucial to understand how glaciers are changing in response to climate change. To date, only few works assess detailed surface elevation change, due to lack of accurate data. In this regard, the near global coverage of multi‐temporal ASTER DEMs (Digital Elevation Model) has a high potential, which is not yet fully exploited. However, ASTER DEMs suffer from low accuracy and significant outliers.

Therefore, novel quantitative methodologies need to be exploited in order to compensate such problems.

In this study, a robust modelling scheme, R‐WLR‐S (RANSAC – iteratively reWeighted Linear Regression ‐ Spline) is developed. Surface elevation change rates and volume change rate are assessed using ASTER DEMs and the SRTM (Shuttle Radar Topography Mission) DEM for two case studies, one glacier in Bhutan and a region at Mt. Cook in New Zealand. Four individual glaciers (Fox Glacier, Franz Josef Glacier, Tasman Glacier, and Murchison Glacier) around Mt. Cook are examined. 4‐dimensional (X, Y, Z, and Time) DEMs are generated by R‐

WLR‐S. Overall, results showed that the glaciers in the Mt. Cook region considered in this study were shrinking, whereas the glacier in Bhutan was stable in general. Moreover, the four individual glaciers monitored in the Mt. Cook region showed that the glaciers on the west side were thickening at the terminuses between 2000 and 2008, whereas started shrinking at a high rate since then.

This study compared R‐WLR‐S with previous methods. R‐WLR‐S could measure the evolution of glacier volume changes in detail, showing a considerable improvement in filtering outliers in original datasets, and fitting elevation change trends robustly.

Front page image: glaciers around Mt. Cook New Zealand, in an ASTER image scene.

(6)

Acknowledgement

I would like to express my great appreciation to all people who helped me along the way.

The master study finally comes to the end, and I am so grateful for everything I have experienced during past two years.

Particularly, I wish to express my gratitude to my supervisor Prof. Andreas Kääb, for his motivated guidance, invaluable suggestions, and especially proposing the topic for me. Ever since I started to work on this topic, it has been a wonderland. I always feel how lucky I am to study at UiO, and being supervised by Prof Andy. Thank you!

Special thanks go to Nagy Robert, hey dude, I made it! We almost kept each other company during the study in first year. Thank you so much my friend, for your help not only in study, but also in life. I enjoyed so much about our trip to China, I hope someday you would visit there again, remember you promised!

I would like to thank my friends and classmates I met in Como as well. Arinze Hawkins, Muhammad Sulaiman, Nigus Kidane, Tran Viet Hoang, I still remember my last night in Italy, we had dinner together in chicken hut! Thank you for helping a rookie who first time living in a foreign country settling down quickly. Good luck with your Phd in USA, Hoang!

I also want to thank Dr. Shangguan Donghui, for accepting me in Chinese Academy of Science. The short three months intern for sure facilitated the development of this thesis, also broadened my knowledge in glaciology. I truly believe I will benefit a lot from this experience!

Heart of gratitude, I am grateful to my parents for never‐ending care and support. For the last six years, so pity that I rarely have chances to go back home and company you. Your love is great and selfless. I cannot express even in words.

(7)

Content

1. INTRODUCTION ... ‐ 1 ‐

1.1 Overview ... ‐ 1 ‐

1.2 Literature Summary ... ‐ 3 ‐

1.2.1 Glacier Elevation Changes from DEM Differencing ... ‐ 3 ‐

1.2.2 Glacier Elevation Changes from Altimetry ... ‐ 5 ‐

1.3 Aim ... ‐ 6 ‐

1.4 Structure of Thesis ... ‐ 6 ‐

2. THEORY ... ‐ 7 ‐

2.1 DEM Generation ... ‐ 7 ‐

2.1.1 Stereoscopic DEMs ... ‐ 7 ‐

2.1.2 Radargrammetric and InSAR DEM ... ‐ 11 ‐

2.1.3 Other Sources ... ‐ 14 ‐

2.2 Glacier Surface Dynamic ... ‐ 15 ‐

3. STUDY AREA AND DATA ... ‐ 17 ‐

3.1 Study Area ... ‐ 17 ‐

3.1.1 Bhutan Himalaya ... ‐ 17 ‐

3.1.2 New Zealand ... ‐ 18 ‐

3.2 Data ... ‐ 20 ‐

3.2.1 ASTER Silcast DEM ... ‐ 20 ‐

3.2.2 SRTM ... ‐ 22 ‐

3.2.3 Glacier Outlines ... ‐ 23 ‐

3.3 Processing Platforms ... ‐ 24 ‐

4. METHODOLOGY ... ‐ 25 ‐

4.1 DEM Co‐registration ... ‐ 25 ‐

4.1.1 Method ... ‐ 25 ‐

4.2 Estimation of Glacier Elevation Change ... ‐ 27 ‐

4.2.1 Robust Modelling Scheme ... ‐ 30 ‐

4.2.2 Data Filtering ... ‐ 32 ‐

4.2.3 RANSAC ... ‐ 34 ‐

(8)

4.2.4 Spline Function ... ‐ 36 ‐

4.3 Estimation of Glacier Volume Change ... ‐ 36 ‐

4.4 Error Budget and Uncertainty Assessment ... ‐ 37 ‐

4.4.1 Regression Method ... ‐ 38 ‐

4.4.2 4D DEM Method ... ‐ 38 ‐

5. RESULTS ... ‐ 40 ‐

5.1 Bhutan Himalaya Case Study ... ‐ 40 ‐

5.1.1 DEM Co‐registration ... ‐ 40 ‐

5.1.2 4D DEM ... ‐ 45 ‐

5.1.3 Glacier Change Analysis ... ‐ 48 ‐

5.2 New Zealand Case Study ... ‐ 52 ‐

5.2.1 DEM Co‐registration ... ‐ 52 ‐

5.2.2 4D DEM ... ‐ 57 ‐

5.2.3 Glacier Change Analysis ... ‐ 60 ‐

5.2.4 Individual Glaciers ... ‐ 63 ‐

6. DISCUSSION ... ‐ 71 ‐

6.1 DEM co‐registration ... ‐ 71 ‐

6.2 Methods Comparison ... ‐ 71 ‐

6.3 Error Estimation ... ‐ 79 ‐

7. CONCLUSIONS AND PERSPECTIVES ... ‐ 81 ‐

8. REFERENCES ... ‐ 83 ‐

APPENDIX A ... ‐ 88 ‐

APPENDIX B ... ‐ 121 ‐

(9)

List of Figures

Figure 1 Multi‐temporal DEM stack and concept of 4‐dimensional DEM. ... ‐ 2 ‐ Figure 2 Example of glacier elevation change extrapolation based on hypsometry. Source: Arendt (2010). ... ‐ 5 ‐

Figure 3 Collinearity conditions for photogrammetry. S1 and S2 are two image acquisition points.

(X1, Y1, Z1) and (X2, Y2, Z2) are space coordinates of two image acquisition points. A represents the same terrain feature. (x1, y1) and (x2, y2) are image coordinated of two image points a1 and a2. ‐ 8 ‐ Figure 4 Stereo imagery generation from nadir looking and backward looking. ... ‐ 9 ‐ Figure 5 Procedure of extracting DEM from ASTER stereo pair using PCI Geomatica... ‐ 10 ‐ Figure 6 The different observation positions and geometry for radargrammetry (Maître, 2010). S1, S2 are the satellites. Bx and Bz are the horizontal and vertical baseline respectively. P represents the ground point. P1 and P2 are seen by two radar images from location S1 and S2. ... ‐ 11 ‐ Figure 7 Interferometric processing of a SAR DEM. ... ‐ 12 ‐ Figure 8 Cross‐track SAR interferometer (flight paths perpendicular into plane). Two SARs fly on parallel tracks and view the terrain from slightly different directions. B stands for baseline, while B stands for effective baseline. From: (Rosen et al., 2000). ... ‐ 13 ‐ Figure 9 SRTM space segments. Shown are the main components of both the X‐ and the C‐band SAR system (Rabus et al., 2003). ... ‐ 14 ‐ Figure 10 Components of the mass balance of a glacier. From: (Cogley et al., 2011). ... ‐ 15 ‐ Figure 11 Ice flux configuration. From: Dr. Andreas Kääb’s lecture note. ... ‐ 16 ‐ Figure 12 Orthoprojected ASTER image of 20 January 2001 showing a section of the Bhutan Himalayas main ridge. The rectangle with letter A indicates the study area corresponding to Figure 13. Coordinates in WGS84 UTM Zone 46N. ... ‐ 17 ‐ Figure 13 Study area corresponding to rectangle A in Figure 12. The 300 m contour lines are derived from SRTM. Glacier extent (red) is generated after manual modification of GLIMS glacier extent database (Armstrong et al., 2005)... ‐ 17 ‐ Figure 14 Histogram of Glacier Bhutan in the year 2000 (based on SRTM). Profile shows a 2697 m drop height between the terminus and headwater. ... ‐ 18 ‐ Figure 15 Orthoprojected ASTER image of 14 February 2002 showing glaciers in Mt. Cook area. The rectangles with letter A indicates the location of Tasman glacier, B indicates Fox Glacier, C indicates Franz Josef Glacier and D stands for Murchision Glacier. The 1000 m contour lines are derived from SRTM. Glacier extent (red) is generated based on manual modification of GLIMS glacier extent database (Armstrong et al., 2005). Coordinates in WGS84 UTM Zone 59S. ... ‐ 19 ‐ Figure 16 Histogram of ‘Glacier Mt. Cook’ study region based on SRTM. Profile shows a 3389 m drop height between the terminus and headwater. ... ‐ 19 ‐

(10)

Figure 17 Glacier Bhutan case study area. Red: GLIMS outlines. Blue: manually corrected, uncompleted glacier coverage due to limitation of ASTER DEM coverage. Background image: ASTER

2006/01/18. ... ‐ 23 ‐

Figure 18 Upper part illustrates the scheme of DEM shift. Lower part shows a trigonometric function between terrain aspect (ψ) and elevation difference (dh) normalized by the slope tangent (tan(α)). The formula in the lower part gives the mathematical relation. The unknown coefficients a, b and c can be solved by a least square fit. From:Nuth and Kääb (2011) ... ‐ 26 ‐

Figure 19 An example of an elevation dependent bias correction using a 3rd order polynomial. From:Nuth and Kääb (2011). ... ‐ 26 ‐

Figure 20 Bias in the SPOT5 DEM due to an anomaly in the satellite attitude and its correction. From:Berthier et al. (2007) ... ‐ 27 ‐

Figure 21 Example of filtering a random pixel in Glacier Bhutan using the 3‐sigma method and extreme elevation change constraints method. ... ‐ 28 ‐

Figure 22 Workflow of the R‐WLR‐S (RANSAC ‐ iteratively reWeighted Linear Regression ‐ Spline) model. ... ‐ 31 ‐

Figure 23 DEM co‐registration step‐by‐step result. Example shows DEM 2001/01/20 co‐registered to SRTM. Description is given in text. ... ‐ 43 ‐

Figure 24 Horizontal shift correction. Left panel indicates an offset, while it has been removed in right panel. ... ‐ 44 ‐

Figure 25 Elevation‐dependent correction (A) and along track correction (B) by fitting polynomial functions. ... ‐ 44 ‐

Figure 26 Histogram comparison before and after each co‐registration step. ... ‐ 45 ‐

Figure 27 Example of a generated 4D DEM in year 2011 for Glacier Bhutan. ... ‐ 45 ‐

Figure 28 Comparison between a DEM generated from the 4D model and an ASTER DEM, both of 2001/01/20. ... ‐ 46 ‐

Figure 29 Hillshade comparison.A: ASTER Image 2001/01/20. B: Hillshade of ASTER DEM 2001/01/20. C: Hillshade of 4D DEM 2001/01/20. Note, only elevations within the glacier outlines shown. ... ‐ 47 ‐

Figure 30 Comparison between the SRTM DEM and a 4D model for the same time. ... ‐ 47 ‐

Figure 31 The original SRTM DEM of Glacier Bhutan with voids. ... ‐ 47 ‐

Figure 32 Hypsometry and rate of glacier elevation change versus elevation in Glacier Bhutan. Red line represents the dh/dt estimated in each elevation bin (100 m). It is plotted with its corresponding error bar (1 standard deviation). ... ‐ 48 ‐

Figure 33 dh/dt map for Glacier Bhutan. ... ‐ 49 ‐

Figure 34 Glacier Thorthormi and Lugge. From: Naito et al. (2012). ... ‐ 50 ‐

Figure 35 Elevation change study for Bhutan site. From: Gardelle et al. (2013). ... ‐ 50 ‐

(11)

Figure 36 Left: ASTER orthoimage anomaly of left branch of studying glacier (2006/01/18).Right:

Landsat 7 SLC‐off image (2006/01/18). ... ‐ 51 ‐

Figure 37 Volume change of Glacier Bhutan between 2000 and 2011 derived from 4D DEM. ... ‐ 52 ‐

Figure 38 DEM co‐registration steps. The example shows DEM 2006/02/09 co‐registered to SRTM. Description is given in the text. ... ‐ 55 ‐

Figure 39 Horizontal shift correction. Left panel indicates an offset, while it has been removed in right panel. ... ‐ 56 ‐

Figure 40 Elevation‐dependent correction (A) and along track correction (B) by fitting polynomial functions. ... ‐ 56 ‐

Figure 41 Histogram comparison for the co‐registration steps. ... ‐ 56 ‐

Figure 42 Example of a generated 4D DEM for the year 2010 for Glacier Mt. Cook. ... ‐ 58 ‐

Figure 43 Hillshade comparison.A: ASTER Image 2012/03/20. B: Hillshade of ASTER DEM 2012/03/20. C: Hillshade of 4D DEM 2012/03/20. ... ‐ 58 ‐

Figure 44 Comparison between a 4D DEM model with an ASTER DEM of 2012/03/20. ... ‐ 58 ‐

Figure 45 Comparison between a 4D DEM with SRTM at the same time. ... ‐ 59 ‐

Figure 46 Original SRTM DEM of Glacier Mt. Cook with voids. ... ‐ 59 ‐

Figure 47 Hypsometry and rate of ice elevation change versus elevation in Mt. Cook. Red line represents the dh/dt estimated in each elevation bin, plotted with its corresponding error bar (1 standard deviation). ... ‐ 60 ‐

Figure 48 dh/dt map of Glacier Mt. Cook. ... ‐ 61 ‐

Figure 49 Error map of regression. The error in each pixel is determined as standard deviation of residuals in the regression. ... ‐ 61 ‐

Figure 50 Volume change of Glacier Mt. Cook between 2000 and 2012 from the 4D DEM. ... ‐ 63 ‐

Figure 51 Individual glaciers. ... ‐ 65 ‐

Figure 52 Overview of elevation change rate for individual glaciers. ... ‐ 65 ‐

Figure 53 (A) Elevation change map of Fox Glacier. (B) Longitudinal profile of elevation and surface‐ elevation thinning rates (green line in A). DEMs are generated from the 4D model and gradually colored. The red‐black gradient line indicates the corresponding dh/dt (right axis). (C)Volume change of Fox Glacier since 2000 from the 4D DEM. ... ‐ 66 ‐

Figure 54 (A) Elevation change map of Franz Josef Glacier. (B) Longitudinal profile of elevation and surface elevation thinning rates (black line in A). DEMs are generated from the 4D model and gradually colored. The red‐black gradient line indicates the corresponding dh/dt (right axis). (C)Volume change of Franz Josef Glacier since 2000 from the 4D DEM. ... ‐ 67 ‐

Figure 55 Slope map of Franz Josef Glacier. ... ‐ 68 ‐

(12)

Figure 56 (A) Elevation change map of Tasman Glacier. (B) Longitudinal profile of elevation and surface‐elevation thinning rates (black line in A). DEMs are generated from the 4D model and gradually colored. The red‐black gradient line indicates the corresponding dh/dt (right axis).

(C)Volume change of Tasman Glacier since 2000 from the 4D DEM. ... ‐ 69 ‐ Figure 57 (A) Elevation change map of Murchison Glacier. (B) Longitudinal profile of elevation and surface‐elevation thinning rates (black line in A). DEMs are generated from the 4D model and gradually colored. The red‐black gradient line indicates the corresponding dh/dt (right axis). (C) Volume change of Murchison Glacier since 2000 from the 4D DEM. ... ‐ 70 ‐ Figure 58 Elevation and dh/dt for selected pixels over Glacier Mt.Cook. The numbered circles in dh/dt map indicate the location of corresponding graphs. The left most point in each graph represents the SRTM. Big blue boxes represent original elevation data in each DEM. Green filled small boxes are points selected by WLR. Red circles are data points selected by R‐WLR‐S. Green line indicates the regression line of WLR, whereas red line represents the regression line of R‐WL‐S.

Yellow line is fitted using a two‐interval quadratic spline function (for deriving the 4D DEM). ... ‐ 75 ‐ Figure 59 Comparison of R‐WLR‐S (R‐IWLR in the Figure), WLR and DEM subtraction method. ... ‐ 76 ‐ Figure 60 DEM subtraction comparison. A: SRTM ‐ ASTER DEM 2012/03/20; B: SRTM ‐ a generated DEM 2012/03/20 from the 4D model; C: a generated DEM (at same time with SRTM) ‐ a generated DEM 2012/03/20 from the 4D model. Red boxes indicate those pixels have invalid dh/dt (> +10 m/a or < ‐30 m/a). ... ‐ 77 ‐ Figure 61 Model comparison. ... ‐ 78 ‐

(13)

- 1 -

1. Introduction

1.1 Overview

Glaciers are considered as one of the best indicators of climate change (Lemke et al., 2007), and are a major contributor to global sea level rise (Meier et al., 2007, Pfeffer et al., 2008, Nuth et al., 2010). Moreover, glaciers are key water sources in many regions (Pieczonka et al., 2013, Bolch et al., 2012). In mountain areas, glacier‐related environmental hazards are a major threat to local population and ecosystem (Kääb et al., 2005). Therefore, understanding glacier evolution becomes imperative key to understanding climate change in cold regions.

Glaciers on Earth are mostly difficult to access. In‐situ glacier monitoring requires major investments in human and material resources. Air‐borne and space‐borne remote sensing techniques enable wide‐range glacier measurement without the requirement of physical access, which provides a huge scientific potential for glacier‐related research. Almost 10 percent of the world's land masses are currently covered with glaciers, mostly in high latitude locations such as Greenland and Antarctica (USGS, 2014). In order to measure such a large number of glaciers, monitoring from space is essential. Nowadays, large amounts of remote sensing data in archives with adequate resolutions and quality provide a substantial contribution to glacier monitoring. In fact, remote sensing monitoring of glaciers achieved comprehensive and encouraging results. For example, the ESA Glacier_cci (climate change initiative) project has developed, evaluated and summarized algorithms for successfully determining glacial parameters such as area, elevation changes and surface velocity fields.

Variations in these glacial parameters are related to each other according to glacier specific characteristics. Among various parameters, glacier elevation is considered to be pivotal, as glacier elevation change rate (dh/dt) is the key factor for inferring glacier mass balance by means of geodetic methods (Bamber and Rivera, 2007). Precise assessment of thinning or thickening rates is a fundamental issue in glacier mass balance studies. DEM differencing is the common method for calculating elevation change over entire glaciers. Several global DEM datasets such as ASTER DEMs and the SRTM DEM make possible to study a large

(14)

number (DEM su et al., 2 these g acquisit disturba distorte Recentl 2012, W multi‐te compar regressi DEM sta to mod to be ex multi‐te elevatio generat

Figure 1

r of glaciers ubtraction) 2010, Pieczo global data tion princip ance or lac ed by side‐

y, robust re Willis et al., emporal DE red with the

ion method acks (Figure el the evolu xplored to emporal AS on and volu ting a 4‐dim

1 Multi-temp

s. Previous to derive e onka et al.,

often are le. For exam k of visual

‐looking ge

egression fr , 2012b, W EMs takes a e DEM subt ds cannot as e 1) not onl ution of gla fit the actu STER DEMs ume change mensional (‘4

poral DEM s

DEM diffe elevation ch

2013, Gar insufficien mple, DEMs contrast on eometry. U

rom multi‐t Willis et al., advantage f traction me ssess the va ly enable ro acier elevat ual elevation

, this thes e from DEM

4D’) DEM fr

stack and c

- 2 - rencing me hange, volu rdelle et al.

nt for indiv s from optic n the groun sing these temporal D

2012a, Me from stacks ethod. How ariation of g obust regre tion change n change tr is aims to:

Ms; (ii) explo rom multi‐t

concept of 4

ethod used ume change

, 2013, Ber vidual DEM cal stereo (e nd, while ra

data effec EMs has be elkonian et s of DEMs, wever, both glacier eleva ession analy e. Appropria

rend. Using : (i) improv ore the vari

emporal DE

4-dimension

the differe e and mass rthier et al., M analysis d e.g. ASTER) adar data (e ctively rem een explore t al., 2013).

with bette DEM subtr ation chang ysis, but also

ate modellin a combina ve the esti ation of ele EM modellin

nal DEM.

ence of two balance (B , 2007). Ho due to the suffer from e.g. SRTM) mains a cha

ed (Nuimura . Regressio er reliability raction and ge. Multi‐te

o make it p ng function ation of SRT imation of evation cha ng.

o DEMs Berthier owever, eir own m cloud can be allenge.

a et al., n from y when robust mporal possible ns need TM and glacier ange by

(15)

- 3 - 1.2 Literature Summary

A large amount of research has been focusing on quantifying glacier elevation change by remote sensing data. The Glaciers section project within the Climate Change Initiative (Glaciers_cci, www.esa‐glaciers‐cci.org) of the European Space Agency (ESA) is an example, which focuses on quantifying glacier area change, elevation change and velocity field. This project systematically evaluated existing methodologies and algorithms for glacier studies based on remote sensing data (Paul et al., 2013), and it can be regarded as methodology guidance for future glacier research.

Based on the source of data, glacier elevation change studies can be divided into DEM differencing and altimetry.

1.2.1 Glacier Elevation Changes from DEM Differencing

DEM differencing is the main method used for measuring glacier elevation changes, due to the large amount of available DEM data. DEM products are discussed in details in chapter 2 of this thesis. DEM differencing mainly contains two approaches; the traditional DEM differencing method uses subtraction of two DEMs (Berthier et al., 2010, Pieczonka et al., 2013, Gardelle et al., 2013, Berthier et al., 2007), and linear regression from multi‐temporal DEMs provides another approach (Nuimura et al., 2012, Willis et al., 2012b, Willis et al., 2012a, Melkonian et al., 2013). Both methods need pre‐processing of DEM data.

It is obvious that DEM data coming from different acquisition approaches and generation procedures will have different characteristics. Methodologies to integrate data from different sources are vital before implementing DEM differencing. In fact, DEM data have been reported to contain multiple systematic errors (Butler et al., 1998, Erskine et al., 2007, Hoja et al., 2006, Karkee et al., 2008, Kääb, 2005, Moore et al., 1991, Podobnikar, 2007, Yang et al., 2011). Such errors could originate from different reasons, such as sensor jitter, topography effect, and inaccurate sensor orientation. Most of the DEM data available to researchers have been processed beforehand, making it difficult to overhaul them in the analysis. Therefore, statistically DEM processing is essential. As for image matching, DEM

(16)

- 4 -

co‐registration is the most crucial step, enabling DEMs to be aligned. Even a small offset in one DEM could easily lead to a wrong estimation in glacier elevation change.

Typical three‐dimensional matching algorithms (3D hereafter) deal with linear translation, rotation and scale (Gruen and Akca, 2005). These algorithms usually use the least square method to minimize the vertical difference of two 3D objects. Nuth and Kääb (2011) developed a slope‐aspect‐based method to co‐register multi‐temporal DEMs. It has been assessed as one of the best algorithms in terms of DEM co‐registration (Paul et al., 2013), with various studies validating the use of this method (Gardelle et al., 2012b, Melkonian et al., 2013, Paul et al., 2012, Gardelle et al., 2013, Pieczonka et al., 2013, Karimi et al., 2012, Nuth et al., 2013).

After DEM co‐registration, elevation differencing can be used to calculate elevation changes. DEM subtraction directly calculates elevation change by differencing two DEMs.

The DEM subtraction method usually cannot provide full coverage due to the limitation of DEM data. Extrapolation is needed to cover whole glaciers and in order to convert elevation change into volume change. This requires knowledge of the glacier area‐per‐elevation bin, also known as hypsometry (Figure 2). To determine changes in glacier volume, the mean or median elevation change (or changing rate) is derived for each elevation bin, and successively multiplied by hypsometry. Using this method, Pieczonka et al. (2013) revealed that most of the glaciers in Aksu‐Tarim Catchment (Central Tien Shan) show a significant surface lowering in the abla on areas of −18.7 ± 5.6 m for the 1976–2009 periods from 1976 KH‐9 Hexagon, SRTM and 2009 SPOT‐5 (French: Satellite Pour l'Observation de la Terre) datasets. Gardelle et al. (2013) comprehensively studied glacier elevation change in Pamir‐Karakoram‐Himalaya using SRTM and SPOT‐5 derived DEMs, contributing to the understanding of a region‐wide glacier evolution. Berthier et al. (2007) obtained an overall specific mass balance of −0.7 to −0.85 m w.e.a^‐1in the Himachal Pradesh using SRTM and a 2004 SPOT DEM.

(17)

Figure Source DEM d metho linear in Cord shows (2012) the Ne Study Southe betwe 2012a) 1.2.2 G Satellit radiati

e 2 Example e: Arendt (2 differencing odology is d regression dillera Darw

that the av used 15 D epal Himala

of Willis e ern Patago en 2000 a ) using ASTE Glacier Elev

te altimetry on or laser

e of glacier 2010).

g from robu described in

(WLR). Me win Icefield,

veraging ele DEMs to rob

aya. The re et al. (2012 nian Ice (S nd 2012. T ER DEMs an vation Chan

y measures r pulses tra

r elevation c

st linear re n chapter 4 lkonian et a , Chile, from evation cha bustly mode esult was in 2b) utilized SPI) Field.

They also s nd SRTM an ges from A surface ele ansmitted f

- 5 - change ext

gression is 4.2 of this t al. (2013) c m multi‐tem

anging rate el glacier el n line with 156 ASTE Results ind studied Nor nd similar m Altimetry

evations by from and r

rapolation b

based on m thesis) and

alculated e mporal ASTE

is around ‐ levation cha

previous st R DEMs an dicate that

rthern Pata methodology

y recording eceived by

based on h

multi‐tempo commonly levation ch ER DEMs an

‐1.6 ± 0.7 m

ange in the tudies abou nd SRTM to SPI is rap agonian Ice y.

the travel a precisely

hypsometry.

oral DEMs ( employs w anging rate nd SRTM. Th m/a. Nuimu e Khumbu r ut Nepal H o study ice pidly losing efield (Willi

time of mi y located a

(detailed weighted e (dh/dt) he result ura et al.

region in imalaya.

e loss in volume s et al.,

crowave altimeter

(18)

- 6 -

(Chelton et al., 2001). Altimetry data, especially ICESat (Ice, Cloud and land Elevation Satellite, (Schutz et al., 2005)) data, has been used for glacial studies (Bindschadler et al., 1989, Moholdt et al., 2010, Bolch et al., 2013, Neckel et al., 2014, Slobbe et al., 2008, Kääb et al., 2012, Kennett and Eiken, 1997). Two approaches have been explored: (i) repeat‐track method, which calculates elevation differences between ground satellite tracks; (ii) cross‐

over method, which interpolates elevation differences at crossing points of ground tracks (Paul et al., 2013). Altimetry data are sparsely distributed on ground, therefore extrapolation is necessary to obtain glacier volume changes, as for the DEM subtraction method.

1.3 Aim

This project aims to precisely map glacier surface elevation change rates (dh/dt), and modelling the variation of elevation change by generating a 4‐dimensional DEM from multi‐

temporal DEMs. The methodology should be able to cope with different errors associated with DEMs. This project will help to better describe how glacier elevation is evolving in time. Moreover, it is expected to be better able to estimate glacier mass balance to a more reliable level compared to other methods.

1.4 Structure of Thesis

This thesis has eight chapters with two appendixes. It begins with the introduction in chapter 1. Chapter 2 presents the theoretical background behind this thesis, including DEM and a brief introduction of glacier surface dynamics. Chapter 3 gives the study data and areas covered in the project. Chapter 4 demonstrates the implemented methodology. The approach of DEM co‐registration and a robust modelling scheme is introduced. Chapter 5 presents an overview of the outcome of this project, and then literally chapter 6 discusses the result. Conclusions are given in chapter 7 and references are organized in chapter 8.

The thesis also contains two appendixes. Appendix A presents all the datasets used in this research. Appendix B contains the code developed for this project.

(19)

- 7 -

2. Theory

2.1 DEM Generation

DEM is a digital model, and a 3D representation of terrain surface, created from elevation data. It has been widely used in geoscientific applications, such as GIS (Geographic Information System), urban planning, and hazard prevention.

DEM can be generated from various sources, including stereo photogrammetry from aerial surveys or optical satellite imagery, interferometry from radar data, LiDAR (Light Detection And Ranging), and topographic maps. Aerial survey uses sensors that are carried on the aircraft, usually airplane or UAV (unmanned aerial vehicle), providing quick and extensive image data acquisition capability. Space‐born stereoscopy and interferometry utilizes sensors that are carried on satellites, to obtain periodical data of earth surface. The repeating orbit enables obtaining multi‐temporal data. LiDAR has the characteristics of fast, accurate and access directly, providing a remarkable ability to obtain terrain information.

2.1.1 Stereoscopic DEMs

Stereoscopic DEMs are generated using photogrammetric principles. Sensors onboard satellites are used to obtain global scale datasets. Existing satellite sensors mainly use two ways to provide DEM generation ability, along‐track stereo such as ASTER, SPOT5, IKONOS, QUICKBIRD, and cross‐track stereo such as SPOT4. Along‐track means stereo pairs are acquired from the same orbit, whereas cross‐track stereo means stereo pairs are acquired from two different orbits.

Stereoscopic DEM uses rigorous photogrammetric solution, including collinearity and coplanarity condition for the conic perspective of a single image line. However, satellite stereoscopy could be more complicated compare with traditional photogrammetry, due to satellite sensor acquisition strategies (Toutin, 2008). As shown in Figure 3, a stereo pair consists of two images of the same terrain area, with a certain overlapping ratio. The image contains displacement information, forming parallax. Elevation can be obtained from parallax after image matching.

(20)

Figure acquisit acquisit image c The Col point), photogr

Where space co length.

Various backwa (e.g. AL this the Advance

3 Collinear tion points.

tion points.

coordinated linearity eq the groun rammetry is

2 3

a f a y

a f a x









Xs, Ys, Zs a oordinate o Coefficients sensors p rd looking LOS). They u esis, ASTER

ed Spacebo

rity conditio (X₁, Y₁, Z₁)

A represen d of two imag quation rep

d object a s:

) (

3 2 3 1

s s s s

X X

X X a



re object s of terrain po

s in the form provide dif

sensors (e.

use the sim is taken as orne Therm

ons for pho ) and (X₂, Y nts the sam ge points a resents the and the im

( ( ( (

3 2 3 1

Y Y b

















pace coord oint; x and y mula are ele fferent acq .g. ASTER)

ilar principl s an examp mal Emission

- 8 - otogramme Y₂, Z₂) are s me terrain f

1 and a₂. e geometry mage points

( )

3 2 3 1

s s s s

Z c Y







dinate of im y are image ements in im quisition st

and forwar le despite t ple to illust n and Refle

etry. S₁ and space coor feature. (x₁

between th s (Figure 3

) ) ) )

s s s s

Z Z Z Z



mage acquis e coordinate

mage rotati rategies. E rd, nadir an their differe rate stereo ction Radio

d S₂ are tw rdinates of t

, y₁) and (x

he sensors ( 3). Collinea

ition point;

e of image p ion matrix.

Examples in nd backwar ent acquisiti oscopic DEM ometer (AST

wo image two image x₂, y₂) are

(image acq arity equati

(1)

X, Y, Z are point; f is th

nclude nad rd looking s ion geomet M generatio

TER) is an i

uisition ion for

e object he focal

dir and sensors tries. In on. The maging

(21)

instrum near in the th genera indepe looking ratio is backw genera

Figure ASTER as the latitud evalua which Individ et al., 2 ASTER packag and ba

ment onboa nfrared (VN hermal infr

ating a ste endent tele

g, together s 0.6 (Abram

ard telesco ate a stereo

e 4 Stereo im GDEM (cur e outcome de and 83 d ated to be a

is automat dual ASTER

2014).

DEM can b ges (e.g. PC ackward im

ard Terra sa NIR, 15 m re

ared (TIR, ereo, with escope asse r sequentia ms, 2000). N ope has on o pair.

magery gen rrently vers of ASTER degrees so round 12 m tically gene DEM is asse

be generate CI Geomatic

ages (3N a

atellite, con esolution), t

90 m reso a spatial emblies. On

lly constitu Nadir telesc nly one ban

neration from sion 2, avail

mission. It uth latitud m (Team, 20 erated by S essed to ha

ed from a s ca, ENVI, an nd 3B) are

- 9 - nsists of th the shortwa olution). Th resolution ne is for na

uting a ster cope consis nd 3B. 3N

m nadir look able in NAS t covers lan

e (Tachikaw 009). GDEM

LICAST soft ve vertical

stereo pair nd ERDAS).

used. Figur

ree separat ave infrared he VNIR su of 15 m dir looking reo (Figure sts of three

and 3B ba

king and ba SA website) nd surfaces wa et al., 2 is a compil tware (www

accuracy ra

using comm For DEM e re 5 demon

te subsyste d (SWIR, 60 ubsystem m

eters. VNIR and anoth 4). The ba bands, 1, 2 and togethe

ckward look ) is a near‐g

s between 2011). Its v

ation of ind w.silc.co.jp/

anges from

mon remot xtraction, o nstrates a ty

ms, the vis 0 m resoluti makes poss

R consists her is for b ase‐to‐heig 2 and 3N, w

er can be

king.

global DEM 83 degree vertical acc dividual AST /en/produc

10 to 60 m

te sensing s only the VN ypical proce

sible and on), and sible for of two ackward ht (B/H) while the used to

product es north curacy is TER DEM ts.html).

(Schenk

software NIR nadir edure of

(22)

extracti of the D

Figure Geomat The pro are load image a coordin System) Point) is how im the pos stereo p betwee of the c can be parallax of two e

ng DEM usi DEM extract

5 Procedu tica.

ocedure star ded then. G and the gro ates can be ), ground s s a referenc

ages are re sition and o pairs. The le n the imag correlation p

extracted x. Moreove eipipolar im

ing the PCI tion using th

ure of extr

rts with def GCP (Ground

ound by ass e obtained surveys, geo

ce point th elated to ea orientation

eft and righ es is along process and

from stere r, a score c mages.

Geomatica.

he PCI Geom

racting DE

fining physic d Control Po sociating im

by a variet ocoded ima at can be c

ch other. M of the sens ht images th

an identica d reduces th

eo pairs by channel can

- 10 - . Al‐Rousan matica softw

EM from A

cal model a oint) determ mage coordi

ty of mean ages, vecto clearly ident Math mode

sor (equatio hen have a al x axis. Us he possibili y matching n be obtaine

et al. (1997 ware.

ASTER ste

and projecti mines the r nates with s, including ors, GIS and

tified in two l utilizes bu on 1). Epip common or ing epipola ty of incorr g two epip ed, reflectin

7) gave a de

reo pair u

on. Image c relationship ground coo g the GPS (G

d topograp o or more i undle adjust olar images rientation. M

r images in ect matche polar image ng the pixe

etailed desc

using PCI

channel 3N p between t ordinates. G Global Posi

hic maps.

images, spe tment to ca s are re‐pro Matching fe creases the es. Followin es and me l matching

cription

and 3B the raw Ground itioning TP (Tie ecifying alculate ojected eatures e speed g, DEM asuring quality

(23)

- 11 - 2.1.2 Radargrammetric and InSAR DEM

DEM generated by radar technique can come from clinometry, stereoscopy (radargrammetry), interferometry and polarimetry methods (Toutin and Gray, 2000).

Among them, InSAR (Interferometric SAR) and radargrammetry based on SAR (Synthetic Aperture Radar) images are mostly investigated.

Radargrammetry technique utilizes the amplitude information in a radar images stereo pair. The principle is similar with optical photogrammetric method. The processing steps for radargrammetric method involve co‐registration of two subset images, matching between co‐registered images, height calculation and then geocoding. The geometry of radargrammetry is demonstrated in Figure 6. Radargrammetry uses the magnitude (intensity) value, and is less affected than phase by atmospheric heterogeneous (Yu et al., 2010).

Figure 6 The different observation positions and geometry for radargrammetry (Maître, 2010). S1, S2 are the satellites. Bx and Bz are the horizontal and vertical baseline respectively. P represents the ground point. P1 and P2 are seen by two radar images from location S1 and S2.

InSAR uses the phase differences between two radar images, acquired with a small base‐to‐

height ratio (Rosen et al., 2000). Figure 8 demonstrates the geometry of cross‐track InSAR.

Typical processing steps for DEM generation from InSAR are illustrated in Figure 7.

(24)

Figure 7 Single lo in brigh height, earth, t unwrap system

7 Interferom ook comple htness (amp

surface def the phase i pped phase

(Geocoding

metric proce ex (SLC) SAR

plitude) on formation, a

information . Finally th g).

essing of a S R data are a radar imag and water v n is unwrap he height in

- 12 - SAR DEM.

arrays of co ges. A SAR vapor distri pped. Then nformation

mplex num interferog bution. Aft n height inf

is transfe

mbers, repre ram contai er removing formation i

rred to req

esenting var ins informa g the effect is calculate quired coo

riations ation of t of flat d from rdinate

(25)

- 13 -

.

Figure 8 Cross-track SAR interferometer (flight paths perpendicular into plane). Two SARs fly on parallel tracks and view the terrain from slightly different directions. B stands for baseline, while B stands for effective baseline. From: (Rosen et al., 2000).

SRTM DEM is a typical InSAR DEM product. The Shuttle Radar Topography Mission (SRTM) is a joint project of the National Aeronautics and Space Administration (NASA), the National Geospatial‐Intelligence Agency (NGA) and the German and Italian Space Agencies (Yang et al., 2011). It was carried on Space Shuttle Endeavour, and flied from 11 to 22 February in 2000 between latitudes 60°N and 57°S, approximately covered 80% of the Earth’s land surface (Rabus et al., 2003). Two antenna pairs operating in C‐band and X‐band were simultaneously illuminating and recording radar signals onboard the shuttle Endeavour (Figure 9). The SRTM DEM has three different levels, SRTM1, SRTM3 and SRTM30, respectively have the 1, 3, 30 arcsec sample spacing. 1 arcsec equals 30 m in horizontal extent. Thus the SRTM1 and SRTM3 are usually referred to SRTM ’30 m’ and ’90 m’ (Yang et al., 2011) .

(26)

- 14 -

Figure 9 SRTM space segments. Shown are the main components of both the X- and the C-band SAR system (Rabus et al., 2003).

2.1.3 Other Sources

There are various other DEM data sources beside optical stereoscopic DEM and radar products. However, most of them are regional, and are not free for public. Among those, for example, LiDAR (Light Detection and Ranging) technology, or Laser Scanning, nowadays is widely implemented in DEM acquisition. LiDAR DEM enjoys high accuracy and free of shadow (Liu, 2008), making it increasingly popular in recent years. It is still developing rapidly in terms of sensor and data processing. LiDAR usually acquires a large amount of points, known as point clouds. Usually, it is necessary to convert point clouds to grid DEM (Agarwal et al., 2006).

(27)

2.2 G

Figure Glacier 10) wh the su mass Accum domin This th this pr differe cannot equati

Glacier Sur

e 10 Compo r mass bala hich involve

m of accum to glacier, mulation is t ated by sur hesis does n

oject. Howe encing) of g t directly r

on:

rface Dyn

onents of the nce represe es both gla mulation an

, while ab typically res rface meltin not explain ever, it is im glacier elev

epresent th amic

e mass bala ents the ma cier dynam nd ablation blation is d

sulting from ng for land‐t glacier dyn mportant to

vation chan he actual g

- 15 - ance of a gl ass change o mics and clim

. Accumula defined as m precipitat

terminating namics in d o keep in mi nge only p glacier mas

lacier. From of a glacier.

mate condi ation is defi s processes

ion (e.g. wi g glaciers.

etail, as thi nd that geo rovides a q s balance.

m: (Cogley e . It is a com

tions. Glaci ned as all s that red

nter snow)

s will be be odetic meas

quantitative This is due

et al., 2011) plex system ier mass ba processes t duce glacie . Ablation i

eyond the s surement (e e interpret e to the co

).

m (Figure alance is that add er mass.

s mainly

scope of e.g. DEM ation. It ontinuity

(2)

(28)

where cross se can res redistrib

Figure 1 A typica glacier mass. C represe represe

is the e ection of th

ult not onl butes mass

11 Ice flux c al example i elevation c Consequentl ent glacier d ent volume c

levation ch he glacier. T

y from acc throughou

configuratio is the redist changes at

ly, unlike ph dynamics. H changes, w

ange rate, b This equatio cumulation

t the glacie

n. From: Dr tribution of a certain re hysical glaci However, if hich in turn

- 16 - b is the ma on shows t and ablatio r system (A

r. Andreas K f glacier ma

egion witho ier models, f integrated n are a reliab

ss balance, hat change on, but also Arendt, 2010

Kääb’s lectu ss, resulting out changin geodetic m d of an enti

ble proxy of

and q is th s in glacier o from the 0).

ure note.

g from ice f ng the amo measuremen re glacier, f glacier ma

he ice flux p r surface ele e flow of ic

flux that can ount of volu nts cannot d elevation c ass changes

per unit evation e, as it

n cause ume or directly changes

.

(29)

3. St

3.1 3.1.1

Figure the Bh area c

Figure lines a modific

tudy Are

Study Ar Bhutan Hi

e 12 Orthop hutan Hima orrespondin

e 13 Study a are derived cation of GL

ea and D

rea imalaya

projected AS layas main ng to Figure

area corresp d from SR LIMS glacie

Data

STER imag n ridge. The e 13. Coord

ponding to r RTM. Glacie er extent da

- 17 - ge of 20 Jan

e rectangle dinates in W

rectangle A er extent (

tabase (Arm

nuary 2001 with letter WGS84 UTM

A in Figure 1 (red) is ge mstrong et a

1 showing a A indicates M Zone 46N

12. The 300 nerated aft al., 2005).

a section of

s the study N.

0 m contour fter manual f y

r l

(30)

- 18 -

The Bhutan study area covers 125 km²area, and is located in the northernmost section of the Bhutan Himalayan main ridge (known as ‘Lunana’, (Gansser, 1970)) with altitudes ranging from 4300 m to 4500 m asl between the main Himalayan divide (Iwata et al., 2004).

Its highest peak reaches to 7300 m (Figure 14). The northern part towards to Tibet Plateau has an approximate elevation range of 5000 m asl with relatively flat terrain. The glacier delineated (‘Glacier Bhutan’ hereafter) has an elevation range from 4329 m to 7036 m according to SRTM DEM. Figure 14 illustrates the elevation profile of the Bhutan study site based on SRTM.

Figure 14 Histogram of Glacier Bhutan in the year 2000 (based on SRTM). Profile shows a 2697 m drop height between the terminus and headwater.

3.1.2 New Zealand

The Tasman Glacier, together with Fox Glacier, Franz Josef Glacier, and Murchison Glacier, are major glaciers around Mountain Cook region in New Zealand (To be called ‘Glacier Mt.

Cook’ hereafter for convenience). Among them, Fox glacier and Franz Josef glacier are located in the west side of Mt. Cook (Figure 15). Tasman glacier is the largest compound valley glacier in New Zealand, which has 28 km length of its trunk and 1 to 2 km in width. It covers an area of approximately 55 km² (220 km² if tributary glaciers are included) (Hochstein et al., 1995). Glaciers in New Zealand are well known as maritime type glaciers (Chinn et al., 2014), affected by high precipitation. Glaciers in this region have rapid response times to climatic perturbations compared with other glaciers in the world.