Brief communication : Rapid machine learning-based extraction and measurement of ice wedge polygons in airborne lidar data

We present a workflow for rapid delineation and microtopographic characterization of ice wedge polygons within high-resolution digital elevation models. The workflow, which is extensible to other forms of remotely sensed imagery, 10 incorporates a convolutional neural network to detect pixels representing troughs. A watershed transformation is then used to segment imagery into discrete polygons. Regions of non-polygonal terrain are partitioned out using a simple post-processing procedure. Results from training and validation sites at Barrow and Prudhoe Bay, Alaska demonstrate robust performance in diverse tundra landscapes. The methodology permits fast, spatially extensive measurements of polygonal microtopography and trough network geometry. 15


Introduction and Background
The objective of this research is to develop and report on a workflow for rapid delineation and microtopographic analysis of ice wedge polygons in airborne lidar data.Ice wedge polygons are the surface expression of ice wedges, a form of ground ice nearly ubiquitous to coastal tundra environments in North America and Eurasia (Leffingwell, 1915;Lachenbruch, 1962).Several thousand ice wedge polygons of typical size may occupy a single square kilometer of terrain, motivating our development of an automated approach to map polygonal boundaries.The key innovation in our method is the use of a convolutional neural network (CNN), a variety of machine learning algorithm, to identify troughs (i.e., polygonal boundaries) within lidar-derived imagery.The same techniques could be applied to any form of remotely sensed data with sufficient spatial resolution; however, using a digital elevation model (DEM) as input, our approach permits the extraction of microtopographic attributes from entire populations of ice wedge polygons at the kilometer scale or greater.
High-resolution inventories of polygon geomorphology are of hydrologic and ecologic interest because decimeterscale variability in polygon microtopography can drive pronounced changes to soil drainage (Liljedahl et al., 2016) and surface emissions of CO2 and CH4 (Lara et al., 2015;Wainwright et al., 2015).Several previous geospatial analyses of polygon geomorphology have been motivated by a recent acceleration in the rates at which low-centered polygons (LCPs) are converted into high-centered polygons (HCPs), a process driven by permafrost degradation which improves soil aeration and stimulates The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-167Manuscript under review for journal The Cryosphere Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.enhanced emissions of CO2 (Jorgenson et al., 2006;Raynolds et al., 2014;Jorgenson et al., 2015;Liljedahl et al., 2016).In these prior studies, the most common form of imagery used to demonstrate this transformation has been historic aerial photography.Thus far, analysis of available datasets has demonstrated a consistent uptick, beginning around 1989, in landscape-scale HCP formation throughout the high Arctic; however, precise rates of geomorphic change have been difficult to quantify, as the surveys primarily have relied on proxy indicators, such as the presence of water in deepening HCP troughs (Steedman et al., 2017).In a related effort to characterize contemporary polygon microtopography, a landcover map of LCP and HCP occurrence across the Arctic coastal plain of northern Alaska was recently developed using multispectral imagery from the Landsat8 satellite at 30 m resolution (Lara et al., 2018).This dataset offers a static estimate of variation in polygonal form over unprecedented spatial scales; however, geomorphology was inferred from the characteristics of pixels larger than typical polygons, preventing inspection of individual features.
Higher-resolution approaches to explicitly segment imagery into discrete polygons have typically been motivated by efforts to analyze trough network geometry.On both Earth and Mars, for example, paleo-environmental conditions in remnant polygonal landscapes have been inferred by comparing parameters such as trough spacing and orientation with systems in modern periglacial terrain (e.g., Levy et al., 2009;Ulrich et al., 2011;Ewertowski et al., 2017).An early semi-automated approach to delineating Martian polygons from satellite imagery was developed by Pina et al. (2006), who first used morphological image processing operations to emphasize polygonal troughs, then applied a watershed transformation (discussed in Section 4.1.3)to identify discrete polygons.This workflow was later applied to lidar-derived imagery from a landscape outside Barrow, AK by Wainwright et al. (2015), but in their experience and our own, robust results at spatial scales approaching a square kilometer or greater were elusive.
Our application of CNNs to the task of identifying polygonal troughs was inspired by the remarkable solutions that CNNs recently have permitted to previously intractable image processing problems.As the performance of graphics processing units (GPUs) has steadily improved over the last decade, CNNs have demonstrated unprecedented skill at tasks analogous to trough delineation, such as cell membrane identification in biomedical images (Ciresan et al., 2012) or road extraction from satellite imagery (Kestur et al., 2018).In training a CNN to analyze lidar imagery from polygonal terrain, our goal was to improve the speed and accuracy of automated polygon delineation in modern tundra environments.Because our method extracts polygons from a high resolution DEM, in enables thorough, landscape-scale inventories of polygonal microtopography, and we anticipate that it will permit precise monitoring of contemporary surface deformation in landscapes covered by repeat airborne surveys.

Study areas and data acquisition
To demonstrate the flexibility of our approach, we applied it simultaneously at study sites near Barrow and Prudhoe Bay, Alaska, two setting with highly divergent ice wedge polygon geomorphology, ~300 km distant from one another.

Barrow
The first study site is located within 10 km of the Beaufort Sea coast (Fig S1A) in the Barrow Environmental Observatory, operated by the National Environmental Observatory Network (NEON).Mean elevation is less than 5 m above sea level, and vegetation consists of uniformly low-growing grasses and sedges.Mesoscale topography is mostly flat but marked by depressions up to 2 m deep associated with draws and drained lake beds (Fig S1B).Microtopography at the site reflects nearly ubiquitous ice wedge development, which becomes occluded in some of the depressions.Ice wedge polygons are of complex geometry and highly variable area, ranging from ~10 m 2 to >2000 m 2 .An airborne lidar survey was flown in August 2012 as part of the U.S. Department of Energy's Next Generation Ecosystems Experiment-Arctic program (https://ngee-arctic.ornl.gov/).The resulting point cloud was processed into a 50 cm horizontal resolution digital elevation model (DEM) with an estimated vertical accuracy of 0.145 m (Wilson et al., 2012).

Prudhoe Bay
The second site is approximately 300 km east of the first and farther inland, located ~40 km south of Prudhoe Bay, AK (Fig S1A).As at Barrow, vegetation consists almost exclusively of low and even-growing grasses and sedges.Mesoscale topography is extremely flat, with a slight (<4%) dip toward the northwest (Fig S1C).Ice wedge troughs are apparent throughout the training site, demarcating polygons of more consistent area than those of Barrow, ~400-800 m 2 .Airborne lidar data was acquired in August 2012 by the Bureau of Economic Geology at the University of Texas at Austin and subsequently processed into a 50 cm resolution DEM (Paine et al., 2015).Vertical accuracy was estimated at 0.10 m.

Polygon delineation algorithm
A high-level summary of our workflow is presented in

Pre-processing
In the pre-processing stage, regional topographic trends were estimated by processing the DEM with a 2D filter, which assigned to each pixel the mean elevation within a 20 m (40 pixel) radius.Microtopography was estimated by subtracting this regional topography from the DEM (Fig S2).In preparation for passing the data to the CNN, microtopography was then converted to 8-bit gray scale imagery.The minimum intensity (0) was assigned to depressions of 0.7 m or greater, and the maximum intensity (255) was assigned to ridges of 0.7 m or greater.

Convolutional neural network
The function of the CNN in our workflow was to identify pixels likely to represent troughs.Conceptually, a CNN is a classification tool which accepts images of a fixed size (in our case, 27×27 grayscale arrays) as input and generates categorical labels as output.The CNN determines decision criteria through training with a set of manually labeled images.The architecture of a CNN consists of a user-defined sequence of components, or layers, which take inspiration from the neural connections of the visual cortex.Our CNN was purposefully constructed with an architecture of minimal complexity, to maximize the efficiency of training and application.Here we briefly summarize the function of each layer in our CNN; for more detailed description, the reader is directed to Ciresan et al. (2012).
The most important components of our CNN were a single convolutional layer, a max-pooling layer, and two fully connected layers (Table 1).In the convolutional layer, a set of 2D filters was applied to the input image, generating intermediate images in which features including concavities, convexities, or linear edges were detected.The max-pooling layer downsized these intermediate images from 27×27 to 9×9 arrays by selecting the highest intensity pixel in a moving 3×3 window with a stride of 3 pixels.Each pixel in the downsized intermediate images was then passed as an input signal to the fully connected layers, which functioned identically to standard neural networks.Two additional components of our CNN were Rectified Linear Unit (ReLU) layers, which enhance non-linearity by reassigning a value of zero to any negative signals output by a preceding layer, and a softmax layer, which converted the output from the final fully connected layer into a probability for each categorical label (i.e., trough or not trough).
During training, the weights of the 2D filters in the convolutional layer and the activation functions of the neurons in the fully connected layers were optimized to correctly predict the labels in a training deck of images.Our workflow was designed to generate the training deck by processing square tiles of manually-labeled imagery with 100 m (200 pixel) edges.
In each of these tiles, trough pixels were delineated by hand in a standard raster graphics editor, a process that required ~1 hour per tile (Fig S3A).Our algorithm imported these tiles, identified the geographic coordinates of each pixel identified as a trough, then extracted the 27×27 neighborhood surrounding each trough pixel from the 8-bit microtopographic imagery.This procedure generated several thousand thumbnail-sized images centered on a trough from each manually delineated tile.
Subsequently, an equal number of pixels not labeled as troughs were selected at random, and the image extraction procedure was repeated.Just prior to training, 25% of the training deck assembled by these methods was set aside to be used for validation.

Polygon extraction
After applying the CNN to classify all pixels at a site as trough or not trough, we extracted discrete ice wedge polygons through application of several standard image processing operations.The first step was elimination of "salt and pepper" noise in the binary image, which we accomplished by eliminating all contiguous sets of trough-identified pixels with an area < 20 m 2 .
Next, we assigned to every non-trough pixel a negative intensity equal to its Euclidean distance from the nearest trough.This created an intermediate image in which each ice wedge polygon appeared as a valley, surrounded on all or most sides by ridges representing the trough network (Fig 1C).At this stage, to prevent over-segmentation, valleys with maximum depths of three pixel-widths or fewer were then identified and merged with the closest neighbor through morphological reconstruction (Soille, 1992).The effect of this procedure was to ensure that the algorithm would not delineate any polygon whose center did not include a point greater than 1.5 m from the troughs.Next, watershed segmentation was then applied to divide the valleys into discrete polygons.Our use of this operation was inspired by its incorporation in the polygon delineation method developed by Pina et al. (2006).Conceptually, this procedure was analogous to identifying the up-gradient region or area of attraction surrounding each local minimum.

Post-processing
In the final stage of delineation, we partitioned out regions of a survey area that had been segmented using techniques described above but were unlikely to represent true ice wedge polygons.For example, polygons were eliminated from the draw in the southern half of the Barrow training site (Fig S1B ), where microtopography is too occluded to permit accurate delineation.Toward this aim, our algorithm tabulated within the boundaries of each delineated polygon (black lines in Fig 1D ) the number of pixels that had been identified by the CNN as troughs (white pixels in Fig 1B).It then eliminated contiguous regions greater than 10,000 m 2 (representing the area of several large polygons) where less than half the boundary pixels of the contained polygons had been labeled as troughs.This procedure had the strengths of being conceptually simple and providing a deterministic means of partitioning non-polygonal terrain from the rest of a survey area.

Microtopographic analysis
To demonstrate our workflow's capabilities for microtopographic analysis, we developed a simple method for measuring the relative elevation at the center each delineated polygon.This served as a rough proxy for LCP or HCP form.In each polygon, we first calculated the Euclidean distance from the trough network of all interior pixels, then identified the pixel most distal from a trough.Elevation was extracted from this pixel and compared with mean elevation in the troughs surrounding the polygon.We thus assigned a number to each polygon which consistently was lower in LCPs and higher in HCPs.

Case study experimental design
We conducted our case study in two phases: calibration and validation.Calibration was performed at two 1 km 2 sites, sampled from the Barrow and Prudhoe Bay DEMs.Leveraging the rapid training and application times of our CNN, we took an iterative approach to extracting polygons from both sites: after training and applying an initial model, new training data were introduced from regions of poor performance to improve skill.After several iterations, our final CNN was trained using manually labeled data from three 100 m tiles at the Barrow training site and one at the Prudhoe Bay site, representing 3% and 1% of the training plots, respectively.The non-trough training data were also supplemented with some additional examples of LCP centers from Barrow, to discourage inaccurate labeling as troughs (Fig S3B).After training the CNN, we tested its extensibility by applying it at two validation sites, also 1 km 2 in area, selected from the same landscapes.Once delineation was complete, to illustrate capabilities for microtopographic analysis, we calculated the relative elevations of polygon centers at the Prudhoe Bay training site, which was chosen for the tendency of its polygons to cluster into areas with high concentrations of either LCPs or HCPs.

Results and discussion
Results from the Barrow training site demonstrate the capability of our workflow to extract polygons of highly variable geometry (Fig 2A).The results also demonstrate success of our post-processing workflow in partitioning out a streambed traversing the southern half of the site and a lakebed in the northeast quadrant.The final delineation makes immediately apparent the complex distribution of polygon type at the kilometer scale.For example, just above the draw in the southwest quadrant is a cluster of comparatively small HCPs that stands in stark contrast to the larger LCPs and intermediate polygons dispersed elsewhere.The complexity of the delineated network underscores the value of an automated approach to polygon extraction, as the manual processing necessary to achieve similar results would require prohibitively large amounts of time.Moreover, the site contains several areas in which larger polygons typically are subdivided by secondary and tertiary ice wedges.Although polygonal boundaries are somewhat ambiguous in such situations, our algorithm provides a consistent, deterministic method of interpretation, which is set by the user's delineation of troughs in the training dataset.
Results from the Prudhoe Bay training site, using the same CNN, demonstrate the ease with which our workflow can The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-167Manuscript under review for journal The Cryosphere Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.be extended across diverse landscapes.We note that, while we had initially trained our CNN using three tiles of manually labeled data from Barrow, successful extension to Prudhoe Bay required incorporation of only one additional tile, from a region in the southeast quadrant characterized by LCPs with exceptionally narrow troughs.The final delineation (Fig 2B ) emphasizes the far more consistent size and geometry of polygons relative to Barrow, which may be associated with the simpler mesoscale topography (Lachenbruch, 1962).The ice wedge network at Prudhoe Bay appears to be less riddled with secondary ice wedges than at Barrow, but some examples are evident in HCPs.The capacity of our workflow for high-resolution microtopographic analysis is demonstrated in estimates of relief at the polygon centers at the Prudhoe Bay training site (Fig 3A).Through relatively simple calculations, we highlight the presence of a distinct cluster of LCPs in the southeast quadrant, which emerges from the mixture of HCPs and intermediate polygons present through the remainder of the site.Using this approach, our procedure provides a means to estimate polygonal geomorphology with unprecedented detail at the kilometer scale, which may improve predictions of soil wetness and polygonscale fluxes of CO2 and CH4 (Lara et al., 2015;Wainwright et al., 2015).Simultaneously, the delineation performed at the same site provides a means to analyze trough geometric parameters, such as spacing and orientation.We acknowledge that, because the surface expression of ice wedges is sometimes subtle or non-existent, ground-based methods remain the highestaccuracy approach to mapping ice wedge networks (Lousada, 2018).Nonetheless, by segmenting machine-delineated polygonal boundaries into individual troughs (Fig 3B), very large datasets can rapidly be generated to estimate statistics at spatial scales unattainable through on-site surveying.

Conclusions
Convolutional neural networks enable rapid delineation of ice wedge networks at unprecedented spatial scales using airborne lidar datasets.Successful application of CNNs is facilitated by employing relatively simple neural architecture, which permits rapid training, testing, and incorporation of new data to improve skill.A single CNN is capable of extracting polygons of highly variable area and geometry from diverse tundra landscapes, using a training workflow that can be completed in a single afternoon.Potential applications for the technology include: generation of high-resolution maps of land cover by polygon type, precise quantification of microtopographic deformation in areas covered by repeat airborne surveys, and rapid Fig 1, using a sample of the data from Barrow.In the first (pre-processing) stage, regional trends were removed from a DEM (Fig S2), generating an image of polygonal microtopography (Fig 1A).Next, the microtopographic information was processed by a CNN, which was trained to use the 27×27 neighborhood surrounding each pixel to assign a label of "trough" or "not trough" (Fig 1B).Each non-trough pixel was then assigned a negative intensity proportional to its Euclidean distance from the closest trough, generating a grayscale image analogous to a DEM of isolated basins, in which the polygonal troughs appear as ridges (Fig 1C).Subsequently, a watershed transform was applied to segment the image into discrete ice wedge polygons (Fig 1D).These steps, and a post-processing algorithm used to remove non-polygonal terrain from the final image, are described in detail below.The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-167Manuscript under review for journal The Cryosphere Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.
The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-167Manuscript under review for journal The Cryosphere Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.Due to the compact architecture of our CNN, training and execution speeds were rapid.Processing information from ~30,000 images, classification accuracies >99% on the training deck and >98% on the validation deck were obtained in less than ten minutes on a single GPU.Application of the trained CNN to a square kilometer of gridded lidar data required ~90 seconds.This speed permitted us to iteratively train the CNN, identify regions of the survey area where it performed poorly, and supplement the training data with labeled imagery from those regions.Optionally, our workflow at this stage also permitted supplementation of the training deck with individual instances of problematic features, such as LCP centers or especially narrow troughs (Fig S3B).This permitted us to increase the skill of the CNN in a targeted manner without the time requirements of manually labeling a full 100 m tile of imagery.
Results from the validation sites demonstrate similar skill in polygon delineation without incorporation of any new training data.At the Barrow validation site (Fig S4A), the algorithm partitioned out another stream bed and three regions where surface inundation prevented accurate delineation of polygonal microtopography.Elsewhere, a mixture of LCPs, HCPs, and intermediate polygons reminiscent of the training site was successfully delineated.At the Prudhoe Bay validation site (Fig S4B), the algorithm partitioned out a thaw lake and the flood plain of a braided stream, in which no ice wedge troughs are present.It also partitioned out a zone just north of the lake where ice wedges may be present in the subsurface, but microtopography is difficult to distinguish.Through the remainder of the area, the algorithm delineated a set of LCPs and HCPs of similar size and geometry to the training site, but with somewhat subtler troughs.
The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-167Manuscript under review for journal The Cryosphere Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.extraction of trough parameters including spacing and orientation.These capabilities can improve understanding of environmental influences on network geometry and facilitate assessments of contemporary landscape evolution in the Arctic.

Figure 1 .
Figure 1.High-level schematic of our workflow.The input to the CNN is 8-bit grayscale imagery representing microtopography (A), estimated by removing regional trends from the lidar DEM (FigS2).The CNN identifies pixels likely to represent

Figure 2 .
Figure 2. Delineation of ice wedge polygons using a single CNN at Barrow (A) and Prudhoe Bay (B) training sites.DEMs are shown on the left, and polygons are randomly colorized on the right.

Figure 3 .
Figure 3. Relative elevations of polygon centers at Prudhoe Bay training site (A).Automated segmentation of the primary trough network at the Prudhoe Bay training site, randomly colorized (B).