LiDAR/DEM development

Airborne laser scanning LiDAR data have a sampling density of 1–2 returns per square meter and a relatively high vertical accuracy (±15 cm). LiDAR has gained widespread use as a tool to generate DEMs that are superior to DEMs from topographic maps or satellite imagery (Nelson and others, Reference Nelson, Reuter and Gessler2009). Previous studies demonstrate that LiDAR data can be successfully used to characterize linear features such as streams and their order, magnitude, width, sinuosity and channel morphology (James and others, Reference James, Watson and Hansen2007; Cavalli and others, Reference Cavalli, Tarolli, Marchi and Dalla Fontana2008; Vianello and others, Reference Vianello, Cavalli and Tarolli2009; Biron and others, Reference Biron, Choné, Buffin-Bélanger, Demers and Olsen2013). Investigations have confirmed that LiDAR data can be used to identify and track glacial crevasses (Arnold and others, Reference Arnold, Rees, Devereux and Amable2006; Jóhannesson and others, Reference Jóhannesson, Björnsson, Pálsson, Sigurðsson and Þorsteinsson2011).

We performed airborne surveys of Haig Glacier on 20 April and 12 September 2015 using a Reigl 580 laser scanner with a dedicated Applanix PosAV 910 Inertial Measurement Unit (IMU). Flight paths were 400–500 m above the glacier surface, at a flight speed of ~50 m s⁻¹ and with a beam divergence of 0.2 mrad, yielding a beam diameter of ~0.1 m. Laser shot densities for the April and September surveys averaged 2.89 and 0.93 m⁻², respectively. We processed LiDAR and flight trajectories using Applanix's PosPac Mobile Mapping Suite which yielded vertical and horizontal positional uncertainties of ±0.15 m (1σ). We generated DEMs used in this study from post-processing point cloud using LAStools (https://rapidlasso.com/lastools/) and the LAStools las2DEM algorithm on classified point clouds. This algorithm triangulates ground classified airborne laser scanning points into a triangulated irregular network from which we created 1 m DEMs using nearest neighbor interpolation. The DEMs were coregistered following the methods described in Nuth and Kääb (Reference Nuth and Kääb2011). This co-registration produced vertical errors of ±0.25 m (1σ) over areas of stable terrain. A full summary of methods pertaining to our LiDAR processing steps and the uncertainties is found in the Supplementary Material of Menounos and others (Reference Menounos2019).

Input layers for crevasse delineation in SOM

Surface crevasses on glaciers form due to brittle failure under tension as ice undergoes extension (velocity divergence) in association with flow over bedrock obstacles, lateral drag from valley walls, changing valley width or orientation, basal sliding, or drawdown of the terminus (Hargitai and Kereszturi, Reference Hargitai and Kereszturi2015). Depending on their genesis, crevasses have different depths, widths and shapes. Most crevasses have nearly vertical walls, and usually propagate downward from the surface. They can widen and deepen under continued extensional stress as well as thermal and mechanical erosion from supraglacial waters that drain into crevasses (Morawski, Reference Morawski2010; Benn and Evans, Reference Benn and Evans2014). The underlying topography and its deformability due to the load of the advancing glacier are both important in their formation zone and evolution (Morawski, Reference Morawski2010), and may cause different internal and external forms of crevasses.

Examination of strain patterns and stress regimes on glaciers through multi-temporal analysis requires a comprehensive description of crevasse morphometry. Mapping crevasses as simple linear features misses some important information that is possible to extract. High-resolution LiDAR DEMs enable the calculation of first, second and third derivatives of local land surface features, which can be used to extract a more complete 3D set of geomorphometric crevasse parameters.

Geomorphometric parameters are topographic variables extracted from elevation models; each can be useful in defining a specific aspect of a landform or shape characteristic. More than 20 geomorphometric parameters have been introduced by different studies (e.g. Florinsky, Reference Florinsky2009; Hengel and Reuter, Reference Hengel and Reuter2008). We examine these parameters based on crevasse shapes in our study area, to select the best factors that delineate crevasses. We consider both cross-sectional (transverse) and longitudinal curvature (Fig. 3), corresponding to the planform and profile curvature of Zevenbergen and Thorne (Reference Zevenbergen and Thorne1987), respectively. Cross-sectional curvature, in the direction perpendicular to the surface slope, defines geomorphological form when surface concavities such as channels and longitudinal crevasses are detected (Wood, Reference Wood1996). Longitudinal curvature is the curvature in the down-slope direction, and helps to identify transverse crevasses. Curvature is important in the morphology of crevasses because of their iconic, relatively sharp cross-sectional profile compared with the rest of the glacier surface (Hengel and Reuter, Reference Hengel and Reuter2008; Hargitai and Kereszturi, Reference Hargitai and Kereszturi2015).

Fig. 3. Schematic illustration of two second-order derivatives of elevation: cross-sectional and longitudinal curvatures. Cross-sectional curvature measures curvature perpendicular to the downslope direction, while longitudinal curvature measures curvature in the downslope direction.

We also identify unsphericity as a factor for further analysis of these features with repeated DEMs. It is a nonnegative curvature, equal to zero on a sphere and positive outside a sphere (Hengel and Reuter, Reference Hengel and Reuter2008); it characterizes how far a surface is from a sphere, which can help in evaluating changes in crevasse rims over time (Fig. 4).

Fig. 4. Changes in the form of rims and troughs of crevasses from aging in the crevasses with no extension or compression dynamics. This evolution can be recognized on DEMs by using the unsphericity parameter. The schematic image indicates the decrease in unsphericity (crisp and sharp edges becoming rounded) over time.

Geomorphometric parameters were examined visually to identify the ones which best exemplify crevasse morphology at our site. They have been displayed and analyzed for crevasse locations to identify parameters that define all parts of a crevasse including its sinuosity, concavity and rim sharpness and consequently its outline and anomalies on the glacier. For extracting these parameters, a local window is passed over the LiDAR DEM and the change at each point relative to surrounding points is derived by a bivariate quadratic function of the surface:

(1)$$Z = \displaystyle{r \over 2}x^2 + \displaystyle{t \over 2}y^2 + sxy + px + qy + f$$

where f is the height at the central point of the window and p, q, r, s and t coefficients are partial derivatives $\left( {p = \; {\textstyle{{\partial z} \over {\partial x}}}, q = \; {\textstyle{{\partial z} \over {\partial y}}},\,\; r = \; {\textstyle{{\partial^2z} \over {\partial x^2}}},\,\; s = {\textstyle{{\partial^2z} \over {\partial x\partial y}}},\,t = \; {\textstyle{{\partial^2z} \over {\partial y^2}}}} \right).$ Parameters greater than the second order do not improve crevasse classification. The optimal kernel size for testing in the methodology has been initially identified visually, selected to give an output layer, which captures the size, orientation and location of individual crevasses.

Crevasses have linear shapes and their footprint compared with the whole study area is relatively low, so we need to use small kernel sizes for better delineation of the elements. However, to reduce high resolution DEM errors that could have significantly affect the curvature parameters of small features (Sofia and others, Reference Sofia, Pirotti and Tarolli2013), through our methodology we have identified 5 × 5 as the best kernel size for all parameters that have been used for crevasses (p, q, r, s and t: quadratic coefficients) as formulated below:

(2)$$ \eqalign{ & {\rm Unsphericity}\,(1/m) = \lsqb {\lcub {r\lsqb {\lpar {1 + q^2} \rpar /\lpar {1 + p^2} \rpar } \rsqb } } ^{1/2} \cr& - t/\lsqb {\lpar {1 + q^2} \rpar /\lpar {1 + p^2} \rpar } \rsqb ^{1/2} \; \rcub \; ^2/\; \lpar {1 + p^2 + q^2} \rpar \cr& + {{\lcub {\,pqr{\lsqb {\lpar {1 + q^2} \rpar /\lpar {1 + p^2} \rpar } \rsqb }^{1/2} -2{\lsqb {\lpar {1 + q^1} \rpar \lpar {1 + p^2} \rpar } \rsqb }^{1/2}s }}}\cr & {{{+\; pqt/{\lsqb {\lpar {1 + q^2} \rpar /\lpar {1 + p^2} \rpar } \rsqb }^{1/2}} \rcub }^2} \rsqb \; ^{1/2}/\lsqb {2{\lpar {1 + p^2 + q^2} \rpar }^{3/2}} \rsqb }$$

(3)$$\eqalign{& {\rm Longitudinal\; Curvature}\; \lpar {1/m} \rpar \cr & \quad = -2\left( {\left( {\displaystyle{r \over 2}p^2 + \displaystyle{t \over 2}q^2 + spq} \right)/\lpar {\,p^2 + q^2\;} \rpar } \right)}$$

(4)$$\eqalign{&{\rm Cross-sectional\;} \,{\rm Curvature}\; \lpar {1/m} \rpar \cr & \quad = -2\left( {\left( {\displaystyle{t \over 2}p^2 + \displaystyle{r \over 2}q^2-spq} \right)/\lpar {\,p^2 + q^2\;} \rpar } \right)}$$

Equation (3) is from Shary (Reference Shary1995), Eqn (4) from Krcho (Reference Krcho1973) and Eqn (5) from Wood (Reference Wood1996). The length of each crevasse is the centerline length, while the width factor is defined from the middle of each feature, perpendicular to the centerline. Orientation, however, is the alignment of the straight line distance between two ends of the feature. We calculate minimum crevasse depth by removing identified crevasses from the DEM and estimating the planar surface through inverse distance weighted (IDW) interpolation using the crevasse edge points. This is a common method in geoscience and is an accurate interpolator for DEM generation with LiDAR data (Su and Bork, Reference Su and Bork2006; Ashraf and others, Reference Ashraf, Hur and Park2017). The lowest point in each crevasse is then found by subtracting the crevasse elevation model and the planar surface that replaces it, giving the deepest detected part of each crevasse in our data (Fig. 5a). These are minimum depths because the LiDAR may not have detected the deepest part of each crevasse, but could be returning an internal reflection from the sides. For accuracy assessment of the results, 50 randomly selected crevasses from different locations of the glacier have been digitized and their extent properties have been extracted manually.

Self-organizing map (SOM) application

ANN algorithms have been widely applied in a variety of different scientific disciplines. Sample applications in geoscience and remote-sensing studies include image classification (Foroutan and others, Reference Foroutan, Kompanizare and Ehsani2013; Maggiori and others, Reference Maggiori, Tarabalka, Charpiat and Alliez2017), crestline mapping (Foroutan and Zimbelman, Reference Foroutan and Zimbelman2017), river flow prediction (Karunanithi and others, Reference Karunanithi, Grenney, Whitley and Bovee1994; Akhtar and others, Reference Akhtar, Corzo, Van Andel and Jonoski2009), hydrological modeling (Dawson and others, Reference Dawson and Wilby2001) and estimation of glacier ice thickness from terrain morphometry (Clarke and others, Reference Clarke, Berthier, Schoof and Jarosch2009; Quiroga and others, Reference Quiroga2013). These networks consist of layers of neurons, which are inter-connected and act as processing units. Each output layer neuron is connected to all neurons in the input layer by synaptic weights (weight vectors). By adjusting the weights of neurons, preferred classification from the network can be obtained.

In general, ANNs consist of layers of neurons which are inter-connected and act as processing units. The SOM or Kohonen network, after its inventor (1982), is one of the most extensively used unsupervised ANN algorithms (Kohonen, Reference Kohonen2001). It executes a characteristic nonlinear projection from high-dimensional space onto a low-dimensional array of neurons, which is useful for identifying and visualizing complex and underlying patterns from small to large datasets. The SOM is a competitive learning algorithm, and competes to give the best representation of the input data during the learning process. SOM has a topology-preserving characteristic (Kohonen, Reference Kohonen2001), meaning that vectors that are close in the input space will be mapped to units that are close in the output map (Vesanto and others, Reference Vesanto, Himberg, Alhoniemi and Parhankangas2000; Bação and others, Reference Bação, Lobo and Painho2005), which is a unique feature. Unlike other ANN algorithms, SOM uses a neighborhood function for identifying the topological properties (geomorphometry here) in the inputs. In addition, unlike with other multivariate methods, with SOM, there is no need for a priori interpolation when data are missing (Richardson and others, Reference Richardson, Risien and Shillington2003).

The SOM architecture classically contains two layers. Apart from the Input layer, SOM has an output layer (map) of nodes (known as neurons) with initial random weights (random images on each of the nodes here). During what is termed the ‘training step’ of SOM, which is an iterative process, input data (geomorphometric layers here) are consecutively presented to the output map. Through this learning procedure, input vectors that have been introduced to the map and the distance from each weight to the sample vector are calculated by running through all the weight vectors. The node in the output map (with its initial random weight) with the minimum Euclidean distance between itself and the input vector is then identified and updated. This node is called the ‘winner’ node or best-matching unit (BMU), and becomes the center. Its neighbors will be updated based on their relative distances to BMU (Fig. 6).The training procedure is formulated in Eqn (5). An input vector is characterized by the vector x = {x ₁, x ₂,…, x _n}, where n is the number of input variables (three geomorphometric parameters here: unsphericity, longitudinal curvature and cross-sectional curvature). The initialized weight vectors connect m × n neurons, so the winner is:

(5)$${\rm Winner} = \arg \lpar {{\min}_{1 \le i \le mn}\; \lcub {\; \parallel w_i(t)-x(t)\! \parallel} \rcub } \rpar \; $$

where ‖ ‖ is the Euclidean distance measure, x(t) is the input vector and w _i(t) is the synaptic weight of neuron connecting output neuron to the input neuron i at iteration t. Selecting the best-input parameters is considered a pre-processing step in this study, as discussed earlier.

Fig. 5. (a) Cross-sectional profile of three crevasses in the LiDAR data (interpolation line in Figure 2). The red line indicates the interpolated surface from an IDW interpolation technique. Subtraction of the measured elevation from the IDW surface for each crevasse gives us the estimated minimum depth of each feature. We use the interpolation surface to filter out the effect of crevasses in calculation of elevation differences on the glacier (i.e. from repeat imagery). (b) Cross-sectional profile of a crevasse in Haig Glacier identified by cross-sectional and longitudinal curvature parameters through ANN methods. (a) Identified as the deepest part of the crevasse, with negative values in both parameters; (b) crevasse walls with opposite signs in both parameters and (c) the rim of the crevasse, with positive signs for both parameters.

Fig. 6. Schematic illustration of the SOM architecture, consisting of input and output layers and weight vectors that connect these two layers. Different gray tones in the SOM layer indicate the degree of the tuning in neighbor neurons, based on their distance to the winner or best matching unit.

Neighborhood is an important concept in SOM, because the weight vector of the winner and output neurons within a neighborhood radius γ of the winner are adjusted in the direction of the input vector as:

(6)$$\omega _i(t + 1) = \omega _i(t) + \alpha (t)\lsqb {x_i(t)\; -\; \omega_i(t)} \rsqb $$

where ω _i(t + 1) is the adjusted weight and α(t) is the learning rate at iteration t. The weights of those neurons outside the neighborhood remain unadjusted. Within the neighborhood nearer neighbors, receive more training, weight adjustment, than farther neurons (Fig. 6). This competitive learning and lateral interaction stage is known as coarse-tuning. During the coarse-tuning stage the learning rate declines gradually from an initial learning rate (α _max) to a final learning rate (α _min), after a total number of iterations (epochs) (t _max):

(7)$$\alpha {\rm (}t{\rm )} = \left( {\displaystyle{{\alpha_{{\rm min}}} \over {\alpha_{{\rm max}}}}} \right)^{1/t_{{\rm max}}}$$

Repeating this process, thus decreases the radius of the neighborhood (γ) gradually during the coarse-tuning stage:

(8)$$\gamma {\rm (}t{\rm )} = \left( {\displaystyle{{\gamma_{{\rm min}}} \over {\gamma_{{\rm max}}}}} \right)^{1/t_{{\rm max}}}$$

It is usually preferred to start with a large initial neighborhood radius, which decreases until the weight of just the winning neuron is attuned in the output layer. This is called the ‘mapping step’. Finally, the output map will become self-structurally shaped by learning, which leads to a topology (regional organizations of neurons) in the output map, and is the fundamental organization of neuron weights that represent separate classes in the input layer. For detecting the existing patterns, a map of the average distance of a neuron to its neighbors can be visualized in the U matrix. This map itself can be organized into different clusters, with a K-means algorithm used in this study. Ultimately, after the SOM is performed, a classified 2-D image is reconstructed from node weights. The value of each pixel in the original map is substituted by the value of its corresponding neuron weights (which is now in a class) and the desired class (crevasses here) on the initial image extent will be detected. SOM structure can be coded in any programing software or conducted in ANN modules of professional software with precise setting options, particularly image processing packages. However, finding the best setting for extracting the target features is challenging, and several settings must be tested and be programed. The SOM in this study was performed in the MATLAB environment.

Considering the linear shape of the crevasses and their small area compared with that of the much larger host glacier, we are technically taking advantage of one of the shortcomings of the SOM algorithm, known as the magnification factor (Cottrell and others, Reference Cottrell, Fort and Pagès1998), to identify linear features such as crevasses. The factor can be defined as the tendency of SOM to overestimate low probability areas, a factor that may not be found in other ANN methods. In this study, the procedure starts on one side, geomorphometry, with selecting the best parameters as inputs. Statistical parameters such as optimum index factor (Chavez and others, Reference Chavez, Berlin and Sowers1982) can be used for this selection. However, in this study, as discussed previously, we carefully chose the best parameters based on our understanding about crevasse geometry. The kernel size for producing the selected parameters from a DEM is considered a part of the SOM configuration in this study. After the parameters are built and selected, each is normalized with a logistic transformation between zero and one.

Assignment of appropriate settings and learning parameters is the most-important step in efficiently extracting a specific feature using SOM in the detail that is needed. The most important control parameters are: (1) the output layer neuron size, which has a direct correlation with the level of details that would be expected for the classification. For the fully delineation of the target feature outline, the optimum size needs to be identified. Initial neighborhood radius also depends on the output layer neuron size. (2) The number of iterations, which is the number of times that the input data will be presented to the map. All these criteria have been selected by trial and error, and the results compared by evaluating the quality of SOM performance using the quantization error (QE) as the criterion:

(9)$${\rm QE} = \displaystyle{1 \over n}\sum \parallel x_i-m_{x_i}\parallel $$

where n is the number of inputs, x _i is the input and $m_{x_i}$ is the BMU. The error describes the neural map fitting to input data, which means Euclidian distance between data vectors and BMU, and the smallest average QE denotes the optimal solution. For example, the necessary number of iterations can be determined by monitoring QE; once it stops decreasing and remains constant, no more are needed. On the other hand, small output layer neuron size may not properly identify the target features and large output size might miss-classify the targets, so for all control parameters, it is critical to identify the optimal and efficient values.

For this study, the random initialization of weights has been used and the initial neighborhood radius adjusted to be equal to the output layer's diagonal size, which is the maximum possible. To identify the optimal map with the best setting and configuration, several sets of configurations have been applied using different configuration parameters. This approach resulted in ~800 classification results, from which repeated and redundant final maps have been removed using principal component analysis as the identification procedure using the methodology introduced by Jolliffe (Reference Jolliffe1972). The final resulting maps have been narrowed down to 20 classified images with the least QE value. Ultimately, the selected maps were vectorized (nonsimplified) and smoothed by Polynomial Approximation with Exponential Kernel in Arcmap before accuracy assessment with a precise manually digitized map. The extent of the manually extracted crevasses, as the reference map, have been examined and compared with the SOM results, in terms of average areas of overlap and outlier. Finally, the most accurate crevasse map has been selected (Fig. 7). In simple terms, the impacts and quality of different setting parameters have been analyzed to gain the best result. Because of the shape of crevasses in this study and precise parameter selection procedure, identifying the target class is straightforward by manual reclassification, but labeling can also be conducted automatically by using unsupervised statistically weight-based labeling procedure (van Heerden and Engelbrecht, Reference van Heerden and Engelbrecht2013). The final crevasse binary raster map has been converted to a vector file and smoothed. Their dimensions were then measured using the final vector image. These procedures can all be coded in a single file, or one can use modules from different software for each step.

Fig. 7. The flowchart indicates the methodology that has been applied in this study to find the best SOM configuration.