Incorporating Aleatoric Uncertainties in Lake Ice Mapping Using RADARSAT–2 SAR Images and CNNs

Abstract

With the increasing availability of SAR imagery in recent years, more research is being conducted using deep learning (DL) for the classification of ice and open water; however, ice and open water classification using conventional DL methods such as convolutional neural networks (CNNs) is not yet accurate enough to replace manual analysis for operational ice chart mapping. Understanding the uncertainties associated with CNN model predictions can help to quantify errors and, therefore, guide efforts on potential enhancements using more–advanced DL models and/or synergistic approaches. This paper evaluates an approach for estimating the aleatoric uncertainty [a measure used to identify the noise inherent in data] of CNN probabilities to map ice and open water with a custom loss function applied to RADARSAT–2 HH and HV observations. The images were acquired during the 2014 ice season of Lake Erie and Lake Ontario, two of the five Laurentian Great Lakes of North America. Operational image analysis charts from the Canadian Ice Service (CIS), which are based on visual interpretation of SAR imagery, are used to provide training and testing labels for the CNN model and to evaluate the accuracy of the model predictions. Bathymetry, as a variable that has an impact on the ice regime of lakes, was also incorporated during model training in supplementary experiments. Adding aleatoric loss and bathymetry information improved the accuracy of mapping water and ice. Results are evaluated quantitatively (accuracy metrics) and qualitatively (visual comparisons). Ice and open water scores were improved in some sections of the lakes by using aleatoric loss and including bathymetry. In Lake Erie, the ice score was improved by ∼2 on average in the shallow near–shore zone as a result of better mapping of dark ice (low backscatter) in the western basin. As for Lake Ontario, the open water score was improved by ∼6 on average in the deepest profundal off–shore zone.

1. Introduction

Lake ice cover is an essential variable for monitoring climate change [1]. Indeed, dates associated with lake ice freeze-up and break-up processes are highly sensitive to changes and variability in air temperature [2]. Lake ice coverage (or ice concentration) also affects weather events and regional climate due to its significant impact on heat and moisture transfer between lake surfaces and the atmosphere [3,4]. Therefore, lake ice coverage products with an accuracy of 10% or better is a requirement identified by the climate community for improving numerical weather prediction and atmospheric reanalysis products [5,6], and for climate monitoring as determined by the Global Climate Observing System (GCOS) [1]. Another operational use needing accurate lake ice cover products is navigation and operations in ice–covered waters. There are many suitable sources of observations available for mapping and monitoring lake ice coverage such as optical satellite data (e.g., Moderate Resolution Imaging Spectroradiometer (MODIS) and Visible Infrared Imaging Radiometer Suite (VIIRS)); however, synthetic aperture radar (SAR) data remain an important source for lake ice mapping at national ice services since they are largely unaffected by the presence of cloud cover and can be acquired regardless of daylight.Using imaging synthetic aperture radar (SAR) for near real–time operational applications, such as sea ice and lake ice mapping, is becoming increasingly possible due to the launch of satellite constellation missions such as Sentinel–1A/B and the RADARSAT Constellation Mission (RCM) in recent years. These missions provide observations at higher temporal resolutions ( daily to weekly coverage depending on latitude) than previously available.

Automated mapping of Great Lakes ice coverage using SAR has been the focus of several studies, beginning with Nghiem and Leshkevich [7]. For the different ice types, the researchers created a library of the co–polarized backscatter signatures at C–band. Knowing the ice conditions and variability in the Great Lakes is of great importance due to their strong impact on the regional climate, ecology, hydropower industry, and navigation [8,9]. In addition to the benefits that SAR observations bring to lake ice mapping, they also fulfill the requirements of relatively high spatial resolution (≤100 m) and large areal coverage (i.e., in ScanSAR mode with nominal swath widths of 250–500 km) of the Great Lakes. Manual production of Great Lakes ice charts by the Canadian Ice Service (CIS) and U.S. National Ice Center (NIC) is operational. This approach is mainly based on visual interpretation of SAR images complemented by auxiliary datasets such as optical satellite imagery and visual observations from ship and aircraft to get insight on historical ice patterns and meteorological conditions [10].

To reduce human biases associated with visual interpretation and to accelerate analysis, several automated/semi–automated ice/water mapping methods have been proposed using SAR. It should be noted that most of these studies have focused on sea ice mapping, which differs subtly from freshwater lake ice mapping. Of all the studies on sea and lake ice mapping from SAR, only a few suggest unsupervised classification, including those by Deng and Clausi [11] and Sobiech and Dierking [12], and Wang et al. [13]. There are also studies that use an unsupervised classification followed by a supervised labeling of ice and water such as those by Leigh et al. [14] and Hoekstra et al. [15].The abovementioned studies have found that the main drawback of only relying on supervised learning is the bias introduced to training sets. This complexity of the scene and emerged variations are caused by the high inter–class and intra–class variability of backscatter [16]. The main advantages of unsupervised learning are noted as not needing the prior analysis/feature extraction of training samples which leads to an easier implementation [12]. Leigh et al. [14] reported an average ice/water classification accuracy of 96.4% using RADARSAT–2 scenes of the Beaufort Sea, Canada and Wang et al. [13] achieved overall accuracy of 90.4% based on RADARSAT–2 scenes of Lake Erie. Quantitative accuracy metric were not given in works by Sobiech and Dierking [12] and Deng and Clausi [11].

Despite relatively limited research on unsupervised classification, there have been a larger number of studies investigating supervised classification for the task of lake ice mapping from SAR. Supervised classification can be applied on backscatter values [17,18]. However, it has been suggested that relying only on the backscatter values results in uncertainties caused by ambiguous correlations between backscatter and ice/water features; for instance, wind roughened areas of water bodies can be misclassified as ice [19]. The addition of other SAR features such as texture has been found to bring better classification accuracy to the model for sea ice mapping [14,20]. To benefit from both unsupervised and supervised learning and mitigate the drawbacks of each approach, Hoekstra et al. [15] used an unsupervised algorithm for region–based segmentation, called Iterative Region Growing using Semantics (IRGS), followed by pixel–based random forest for labelling. In their research, the potential of the operational usage of the proposed automated method for lake ice cover mapping was demonstrated as a high average ice/water classification accuracy of 95.8% was achieved using RADARSAT–2 scenes of Great Bear Lake (Canada).

Another stream of classifying ice pixels using SAR are methods that focus on learning the features directly from the data [21,22,23]. These methods are based on machine vision techniques that emulate human visual system’s ability to identify and localize objects while assimilating information at different spatial scales based on prior knowledge. The most popular framework for this task in image processing is a convolutional neural network (CNN). The work presented herein builds on previous studies by Wang et al. [23] and Scott et al. [21] to address the challenge of retrieving ice and water from SAR as a result of similar backscatter signatures for ice and water in lakes. A CNN architecture is proposed that incorporates aleatoric uncertainty estimates. Algorithm settings such as input patch size and the training dates are inferred from successful experiments in the aforementioned studies. Improvement in accuracy when including bathymetry within the classification algorithm is also evaluated. Finally, potential applicability of the proposed method to other applications and suggestions to further improve the mapping accuracy are discussed.

2. Background

CNN applications via remote sensing imagery have shown promise in recent studies, having been successfully applied to optical imagery and SAR observations for classification, object detection, and data fusion tasks. These studies have implemented CNN as either the core of their deep learning system or a significant part of them. Lee and Kwon [24] have demonstrated how CNN can be used for target detection using hyperspectral images, and Jiang et al. [25] and Kang et al. [26] have applied Gabor–Deep and contextual region–based CNN on SAR images, respectively, to detect objects. To improve CNNs’ performance in fully retrieving the spatial information, new architectures and expansions have been proposed [27]. As an example, capsule networks (CapsNets) as new neural network structures with the capability of recognizing the spatial information of the image space are being merged with CNNs to classify complex image scenes [27].

Few studies, however, have been conducted on the application of CNNs to the problem of lake and sea ice mapping. Wang et al. [23] mapped ice concentration in the Gulf of Saint Lawrence using RADARSAT–2 images. They compared ice concentration using different methods of multi–layer perceptron (MLP40), CNNs, and an ice concentration product derived from passive microwave data using an algorithm developed based upon ARTIST project (Arctic Radiation and Turbulence Interaction STudy) known as ARTIST Sea Ice (ASI) [28]. Boulze et al. [29] used CNNs to classify sea ice types using 2018 and 2020 Sentinel–1 SAR images and gained an average classification accuracy of 91%. CNNs for lake/sea ice classification can be used in a baseline structure such as in the work reported by Scott et al. [21] and Wang et al. [23], or they can be utilized in complex architectures that have shown better performance in specific tasks such as work by Tom et al. [22] that uses a semantic segmentation network proposed by Chen and Jahanshahi [30] for classifying ice/snow–covered and open water areas in alpine lakes of Switzerland.

In the research of Scott et al. [21], which builds on Wang et al. [23], CNNs are used, and quality control measures are provided based on the uncertainty of the CNN output with the goal of creating reliable sources for incorporating lake ice information in operational data assimilation. In these studies, CNNs have been shown to have powerful feature learning capabilities; learning the features that map the image data to a set of training labels that come from sea ice and lake ice charts. Such training labels introduce potential error and uncertainty into the CNN model predictions, which is due to inherent errors in both spatial representativeness and image interpretation of ice experts using SAR imagery and other auxiliary data.

In this paper, image analysis charts are used to train a CNN to predict ice and water in Lake Erie and Lake Ontario (Canada/United States of America). Given the uncertainty in the ice chart labels, a previously proposed framework [31,32] is adopted to both quantify the uncertainty and obtain improved model predictions by taking this uncertainty into account. The important aspect to consider here is that aleatoric uncertainty cannot be reduced by inputting more data or training a deeper network. Aleatoric uncertainty arises from the nature of data; in other words, it is a measure used to identify the noise inherent in data. This uncertainty quantification can provide in–depth insight into errors inherent in the observations and help us refine where more confidence is expected in classification results. Aleatoric uncertainty is in theory different from uncertainties that are associated with model parameters, called epistemic uncertainty, although there is some ambiguity in trying to separate the two [33].

3. Study Area and Datasets

3.1. Study Area

The study area includes Lake Erie and Lake Ontario, two of the Laurentian Great Lakes of North America (Figure 1). Lake Erie covers an area of 25,655 km${}^2$, with an average depth of 19 m. Lake Ontario with a surface area of 18,960 km${}^2$ has an average depth of 86 m. The Great Lakes do not typically completely freeze over during ice seasons, although Lake Erie sometimes does due to its shallow depth. During winter 2014, because of persistent low temperatures, the Great Lakes recorded their second largest ice coverage since 1973 (NOAA/GLERL NOAA Great Lakes Environmental Research Laboratory–Historical Ice Cover [34].

Hourly air temperature measurements from four weather stations adjacent to the Lake Erie shoreline and five weather stations adjacent to Lake Ontario were used to calculate the average temperature over the lakes during the study period (11 January to 4 April 2014). These temperatures, in addition to the ice analysis charts originally used for model training and testing, provided additional support for the model prediction maps. The average hourly air temperature measurements of Lake Erie and Lake Ontario are shown in Figure 2. Air temperatures corresponding to the times within the hour of SAR acquisitions are also shown in this figure.

3.2. SAR Imagery

RADARSAT–2 images (HH and HV polarizations) acquired on 29 dates during the ice season of 2014 were used for this study. All images were obtained in ScanSAR Wide beam mode with a pixel spacing of 50 m and a scene size of 500 km × 500 km in range and azimuth. This mode of image acquisition is used in operational ice monitoring due to wide swath and appropriate spatial resolution for this task. The temporal resolution is 1–3 days and incidence angle ranges from 20${}^{\circ }$ to 49${}^{\circ }$ for both Lake Erie and Lake Ontario. The list of dates and acquisition times as well as information regarding the dataset split into train, test, and validation, can be found in

Table 1. List of RADARSAT–2 image data used to train, test and evaluate the model.

Acquisition Date	Acquisition Time	Lake Erie			Lake Ontario
		Train	Test	Evaluation	Train	Test	Evaluation
11 January 2014	23:15
12 January 2014	11:26
14 January 2014	23:27
15 January 2014	11:40
18 January 2014	23:11
19 January 2014	11:22
22 January 2014	11:35
28 January 2014	23:19
29 January 2014	11:31
4 February 2014	23:15
5 February 2014	11:26
7 February 2014	23:27
14 February 2014	23:23
21 February 2014	23:19
22 February 2014	11:31
25 February 2014	11:44
28 February 2014	23:15
1 March 2014	11:26
3 March 2014	23:27
4 March 2014	11:40
7 March 2014	23:11
8 March 2014	11:22
10 March 2014	23:23
18 March 2014	11:31
20 March 2014	23:31
25 March 2014	11:26
28 March 2014	11:40
1 April 2014	11:22
4 April 2014	11:34

. Of the 29 images that cover Lake Erie and/or Lake Ontario, ten are used for training and ten are withheld as independent test data. The trained model is evaluated on 19 dates in total (due to scene coverage six are used to evaluate ice mapping in Lake Ontario, and all 19 images are used to evaluate ice mapping in Lake Erie).

RADARSAT–2 images (HH and HV polarizations) acquired on 29 dates during the ice season of 2014 were used for this study. All images were obtained in ScanSAR beam mode with a pixel spacing of 50 m and a scene size of 500 km × 500 km in range and azimuth. This mode of image acquisition is used in operational ice monitoring due to its wide swath and appropriate spatial resolution for this task. The temporal resolution is 1–3 days, and the incidence angle ranges from 20${}^{\circ }$ to 49${}^{\circ }$ for the two lakes. The list of dates and acquisition times as well as information regarding the dataset split into Train, Test, Validation, can be found in. Of the 29 images that cover Lake Erie and/or Lake Ontario, ten are used for training and ten are withheld as independent test data. The trained model is evaluated on 19 dates in total (due to scene coverage six are used to evaluate ice mapping in Lake Ontario, and all 19 images are used to evaluate ice mapping in Lake Erie).

3.3. Image Analysis Charts

Image analysis charts provided by the Canadian Ice Services (CIS) were used to train the model and validate the lake ice product. Image analysis charts are manually generated analyses of SAR imagery in which homogeneous regions are outlined and assigned a label corresponding to the overall fraction of ice concentration in delineated polygons and partial fractions of up to three different ice types. Unlike daily ice charts, image analysis charts are based on the interpretation of SAR images only and have a valid time that corresponds to that of the image acquisition. For labeling water and ice, all pixels with ice concentration less than 10% were labelled as ice. This labeling criterion was done to ensure that labeled pixels were representative of ice and water and other pixels that do not meet the criteria were not included in the training and testing of the CNN. Two examples of image analysis ice charts total ice concentration values (the variable used for training the model) on 19 January and 1 April 2014, are shown in Figure 3. These sample images are taken from the freeze–up and the breakup periods. Freeze–up corresponds to the period of time between the beginning of ice formation and the formation of a complete ice cover (or maximum ice extent in the case when a lake does not form a complete cover) while break–up refers to the period of time between the onset of spring melt and disappearance of ice from the lake surface.

3.4. Data Processing

Due to the coherent nature of the imaging process, SAR images are contaminated with speckle noise. To reduce this noise, a non–overlapping 4 by 4 averaging filter is applied to each image. Each pixel in the resulting image represents the backscatter of a 200 m by 200 m region. This pre–process increases the computational efficiency as well by reducing the data volume. To identify ice and water pixels, the CNN requires image patches. Patch size and stride counts were selected based upon an approach similar to that of Wang et al. [23] to opt the suitable size for such application regarding the original spatial resolution. Patches of 20 by 20 pixels (4 km × 4 km) were extracted with a stride of 3 pixels between the center of each patch. Image analysis charts were interpolated to SAR grids using nearest neighbour interpolation and then same patches were recorded.

4. Methodology

4.1. CNN Classification with a Custom Loss Function for Aleatoric Uncertainty Estimation

A typical CNN may receive as input an image patch and map that patch to a single output. For an example of this type of architecture for lake ice classification see Scott et al. [21]. To learn the uncertainty in the data, the approach taken here is to use a CNN architecture that can predict two items, $y_{predict}^W$ and $y_{var}^W$, with the second one for the uncertainty (represented here as a variance); W represents the network parameters. To do this for our binary classification problem, we followed Kendall and Gal [36] and perturbed the logit produced by the CNN ($y_{predict}^W$) by a sample drawn from a Gaussian distribution with a mean of zero and variance that is a learned parameter ($ϵ\sim N(0,{(y_{var}^W)}^{0.5}$). Each perturbed logit was processed by a sigmoid function ($P_i=sigmoid(y_{predict}^W+ϵ$)) followed by a binary cross–entropy loss function as

$$loss=-{\Sigma }_i(1/T{\Sigma }_t(y_{i}log\left(p_{i,t}\right)+(1-y_i)log\left(p_{i,t}\right))$$

where T is the number of Monte Carlo runs (set as 50 in our algorithm) for Gaussian sampling, $y_i$ is the label value of for the ith sample, and $p_{i,t}$ is the CNN output for the ith sample for the tth Monte Carlo run. The mean value of these binary cross–entropy loss functions was then used to update the network weights and the variance. Note that while the mean value could be compared with the label (supervised learning), the variance was updated through backpropagation (unsupervised learning).

The learned variance was for the logit, but for the end–user a predicted probability of ice and its representative variance at each lat/lon value was more useful as it tells the user directly about the uncertainty in the desired quantity. The CNN model could predict a logit and a variance of that logit at each lat/lon value. To obtain our desired quantities, we generated a set of ice probabilities by perturbing the logit with a perturbation drawn from a Gaussian distribution with the variance equal to that specified at the given lat/lon value ($y_{predict}^W+N(0,y_{var}^W))$. Each perturbed logit was passed through the sigmoid and the mean and variance of these predicted probabilities was computed as predicted probability of ice and its representative variance.

4.2. Incorporation of Bathymetry

Ice forms as lakes begin to cool down. Air temperature, wind speed, lake size, and lake depth drive processes that lead to ice growth [37]. Lake surface area is a variable that relates the time to reach isothermal conditions as the surface heat flux is directly related to the lake surface area. Lake depth also contributes to cooling and freezing time as the deeper a lake is the longer it generally takes to cool [38]. Given the links between freeze–up and bathymetry, we experimented adding bathymetry as an auxiliary data source to boost the CNN performance. To avoid learning unphysical correlations that could be enforced if detailed bathymetry information is included, the bathymetry is classified into different levels. Zonal bathymetry is created based on classification thresholds proposed by Riseng et al. [39]. Considering the depth thresholds that describe each zone, instead of using bathymetry values, we use assigned values of 1 to 4 for Lake Erie and 1 to 5 for Lake Ontario (Figure 4). Riseng et al. [39] identified five bathymetric thresholds for the Great Lakes that are related to processes such as littoral energy, aquatic vegetation extent, stratification limits, and light extinction. These zones are identified as: coastal margin zone (0–3 m depth zone), shallow near–shore (3–5 m), deep near shore (5–15 m for Erie; 5–30 for Lake Ontario), shallow off–shore (15–30 m; Lake Erie only), deep off–shore (30–100m), and profundal off–shore (>100 m; Lake Ontario only). These bathymetry maps coincide with visual patterns we have observed with ice formation in these lakes, in that the ice normally forms first in the western basin of Lake Erie, and near the shoreline in Lake Ontario.

4.3. CNN Architecture and Hyperparameter Selection

The CNN architectures used were based on Scott et al. [21], with some modifications. The first architecture tested, to provide a baseline, was identical to the one suggested by Scott et al. [21]. To incorporate aleatoric uncertainty, this architecture was modified to incorporate two outputs after the final dense layer instead of one (Figure 5a). The sigmoid function that previously followed the dense layer was removed as this was incorporated in the custom loss function. Therefore, the base CNN architecture predicted two quantities, a logit and a variance. Python was used as the programming language. Key utilized libraries include Tensorflow, Tensorflow probability for developing CNNs and NetCDF4 and scipy for handling data processing and post processing steps. To run the codes, computational resources of two v100 GPUs, 18 CPUs with a total of 90 G of memory has been utilized.

To evaluate the impact of including bathymetry, the baseline CNN was first trained and evaluated (without incorporating aleatoric uncertainty). In this case the architecture proposed by Scott et al. [21] was used, except having an additional input, which consisted of zonal bathymetry. Bathymetry is not at the SAR image pixel scale. A zonal bathymetry class is assigned to each lat/lon value on the label grid. Therefore, the classes are incorporated in the CNN at a much lower level when the CNN output is of the same dimension as the label (Figure 5b).

Finally, to incorporate both bathymetry and aleatoric uncertainty the same CNN architecture was used, but with two outputs after the final dense layer, and with the removal of the sigmoid and the use of custom loss function. In this setting however, the way the bathymetry data are interpreted by the CNN is different than without aleatoric uncertainty. These data are part of the Monte Carlo sampling that is used to find the optimal value of the logit and variance, and the weights linking the bathymetry to the network output (mean and variance) are updated in the backpropagation. In this way, the aleatoric uncertainty can be a function of not only the input image patch, but also the bathymetry.

To obtain a smooth reduction in the loss function when both bathymetry and aleatoric uncertainty are incorporated, transfer learning was used. Instead of initializing by random weights, we used weights of trained network with aleatoric loss when no bathymetry was incorporated. The part of the network with transferred weights is delineated by a red box in Figure 5. This method was found to lead to a smooth decrease in the loss function. Before applying this technique, other popular methods for tuning hyperparameters and optimizers for weight initialization and updates were tested but showed to have a minimal impact on the loss trend. We also tested the orthogonality of the weight matrices while training because this can contribute to an uneven decrease in the loss function, we monitored the Forbenius norm of weights (W) calculated by W${}^T$W, and found no sign of losing the orthogonality of the weight matrices while training.

For all CNN architectures, hyperparameter tuning was carried out to ensure an effective learning strategy. We tested different optimizers such as Adam and stochastic gradient descent (SGD) and tuned respective parameters. We found SGD to perform the best for our problem. After selecting SGD, the learning rate was tuned. A cyclic learning rate schedule was found to be optimal in helping the model converge compared to step decay [40].

5. Results

Using typical scores utilized in forecast verification of binary predictands, results are evaluated to infer each model’s performance. Four models are compared: (1) CNN baseline model; (2) CNN with zonal bathymetry; (3) CNN with aleatoric loss; and (4) CNN with aleatoric loss and zonal bathymetry. All models use both the HH and HV SAR images as input data. Contingency tables (

Table 2. Contingency table and prediction scores used for model evaluation.

Model Predictions	Ice Chart
	Ice	Water	Marginal Total
Ice	a	b	a + b
Water	c	d	c + d
Marginal Total	a + c	b + d	a + b + c + d = n
proportion correct ice = a/(a + c) proportion correct water = d/(b + d) total proportion correct = (a + d)/n missed ice = c/(c + d) false alarm (missed water) = b/(a + b)

) for Lake Erie and Lake Ontario for total number of 8,644,181 and 2,885,246 pixel counts are shown in

Table 3. Reported metrics averaged on test dates in Lake Erie.

Model	Proportion Scores
	Correct Ice	Correct Water	Total Correct	Missed Ice	Missed Water
Baseline CNN (inputs: HH-HV)	0.991	0.468	0.967	0.293	0.025
CNN with bathymetry (inputs: HH, HV, bathymetry)	0.992	0.400	0.969	0.179	0.025
CNN with aleatoric loss (inputs: HH, HV)	0.992	0.467	0.968	0.261	0.025
CNN with aleatoric loss and bathymetry (inputs: HH, HV, bathymetry)	0.996	0.421	0.970	0.160	0.027

and

Table 4. Reported metrics averaged on test dates in Lake Ontario.

Model	Proportion Scores
	Correct Ice	Correct Water	Total Correct	Missed Ice	Missed Water
Baseline CNN (inputs: HH, HV)	0.965	0.763	0.819	0.017	0.390
CNN with bathymetry (inputs: HH, HV, bathymetry)	0.910	0.870	0.881	0.038	0.272
CNN with aleatoric loss (inputs: HH, HV)	0.949	0.760	0.812	0.025	0.397
CNN with aleatoric loss and bathymetry (inputs: HH, HV, bathymetry)	0.918	0.840	0.862	0.036	0.313

, respectively. Figure 6 presents results for 15 January 2014, in Lake Erie and Figure 7 shows results of 25 March 2014, in Lake Ontario. Although only one prediction map is shown for each lake, we have observed the same patterns on different dates, and our interpretations are transferable to those dates. Scores in Table 3 and Table 4 are quite close with slight improvements when using bathymetry, and bathymetry and aleatoric uncertainty. However, visual comparisons show substantial improvements in specific regions. For instance, visual comparisons show improvements in mapping ice in the western basin of Lake Erie (Figure 6d,g) when bathymetry is added with/without aleatoric uncertainty, respectively.

In Lake Erie, when bathymetry information is included, scores from scenes with a large water coverage confirmed that the ability of the model to predict open water in the central portion of the lake is significantly reduced leading to a high number of water points misclassified as ice (shown with a solid red line (labelled as region 2)). Dropped performance in central regions is the expense of using large scale zonal bathymetry whereas same values covers different backscatter signatures. The bathymetry and aleatoric uncertainty together provide a good estimate of water points and ice points, with larger aleatoric uncertainties near the southern shore of Lake Erie. This can also be seen in Figure 6, comparing Figure 6c with Figure 6g in the west basin. Results also show that aleatoric uncertainty is high around the ice edges. We also observe that aleatoric uncertainty is generally correlated with small scale image features with less than 2 km by 2 km area. Notice also that there is a patch of smooth dark ice in the western basin shown with a red solid line ellipsoid in Figure 6a (region 1) that is incorrectly classified as water in CNN baseline (Figure 6c). When aleatoric uncertainty is included (Figure 6e), this region is still misclassified, but it is also associated with high uncertainty. When both bathymetry and aleatoric uncertainty are included, the classification is improved and the uncertainty in this region is reduced. At the same time, when both bathymetry and aleatoric uncertainty are included, the patch of smooth ice is identified appropriately as an area associated with high uncertainty. We also observed that models were not successful in identifying ice shown with a red dashed line circle in Figure 6a (region 3). Lack of visible texture in the SAR image is the main reason for classifying ice as water.

In both Lake Erie and Lake Ontario, CNN is not able to fully differentiate open water from ice under wind–roughened water conditions. In Lake Ontario, this effect is more visible due to the larger region of open water (reduced ice cover) and can be seen in the central parts of the lake for the image considered here (Figure 6c). By adding bathymetry, the number of misclassified ice pixels in the central portion of the lake decreased. Using aleatoric loss and bathymetry further decreases these misclassifications. A notable challenge in all applied models is that the ice next to the northern part of the Lake Ontario is consistently misclassified as water (indicated by the red ellipsoid with pointing arrow in Figure 7a). Ice charts and air temperature profiles before and after this date indicate that this region is covered by thin new ice. This region is dark and lacks visible texture in the SAR image, and it is therefore not surprising that it is classified as water. A similar feature can be seen in the western basin of Lake Erie shown with a red solid line ellipsoid in Figure 6a (region 1), as explained earlier. Therefore, the proposed approach is not able to overcome this problem in both Lakes.

All models provided reasonable accuracy similar to other research in this field; however, we noticed nominal differences among different models’ performance. This initial finding motivated us to seek deeper insight by comparing models’ performances within lakes in each zonal bathymetry class. As this classification can better represent lakes’ ecology compared to sub–basin boundaries, zonal bathymetries were chosen as a classification scheme for comparing the scores. Results show that in Lake Erie, using aleatoric loss and including bathymetry in shallow near–shore zone (zone 2) improves the ice score by ∼2 percent in average. This is an interesting finding as it shows the CNN with aleatoric loss and bathymetry can partially overcome the above–mentioned challenges in detecting dark ice in the western basin. In Lake Ontario, in the deepest zone, profundal off–shore (zone 5), using aleatoric loss and including bathymetry improves open water score by ∼6 percent in average. We have also observed that water pixels in shallower zones leads to much lower scores compared to deep zones in both Lake Erie and Lake Ontario. Thus, we can use baseline CNN in deeper zones with more confidence, but for shallow zones (Zones 1 and 2), relying on the baseline CNN decreases the accuracy of mapping water pixels dramatically.

Here, we have shown the predictions for only one date per lake. However, we have observed the same patterns on different dates, and our interpretations are transferable to those dates.

6. Discussion

In this paper, the incorporation of aleatoric uncertainty in a CNN was elaborated and results presented. The model was designed to produce two outputs to estimate the probability of ice and water, and the variance of this estimate. Adapting the model by using the custom loss function provides uncertainty estimates of the predicted probability that are found highly effective in retrievals, specifically in sections of the lakes that exhibit more complex textures such as along lakeshores and ice edges. To gain deeper insight into how uncertainty informs predictions, one possible method would be to examine how misclassification correlates with higher levels of uncertainty. This will frame an approach to improving predictions either by creating a better training set or by fusing the predictions in focused areas with other observations/model predictions. Due to limited access to an independent source of observations with equivalent/higher spatial resolution for the same acquisition dates, such experiment could not be performed in our research. One potential candidate for running the experiment was Landsat images, but images covering the study area were obscured by clouds, making them an inconsistent and thus unsuitable source of information for this study. There are other potential data sources that are discussed hereafter.

Despite the high spatial resolution and imaging capabilities of SAR under cloudy conditions, interpretation of SAR imagery is challenging and requires expert knowledge. Therefore, the preparation of training labels for supervised learning algorithms is costly, as is the case for lake ice classification. For this research, we used available ice analysis charts to provide our training labels. These charts are created at CIS by experts with decades of experience interpreting SAR imagery. Visual interpretations and simultaneous mapping are done by looking at the SAR images as well as historical ice patterns and meteorological conditions. These labels do not contain small–scale features in the ice cover, and the richness of visual cues across ice concentration levels is less than that of sea ice concentration datasets. This finding, in combination with the fact that a deeper CNN was not useful (not shown), leads us to hypothesize that the information content of image analysis charts for lake ice is lower than that for sea ice. We envision creating an enhanced training set using other data sources, such as optical imagery with a high revisit time that increases the chances of obtaining cloud free images that can also be incorporated in the classification framework (uncertainty characterization). A practical example of such a potential data source is the Daily Lake Ice product generated by the Climate Change Initiative (CCI+) Lakes project of the European Space Agency through processing of MODIS Terra/Aqua top–of–atmosphere reflectance data (MOD/MYD02). However, the prospective research is certain to face challenges, such as lower spatial resolution than SAR and cloud obstruction. By understanding the spatio–temporal variability of ice cover, we can confine the use of the lower resolution observations to regions with less variability. On the other hand, more frequent observations from MODIS will allow us to trace pixels of ice/water and fill the gaps caused by cloud cover. However, we would have to quantify any uncertainties introduced by a composite product.

7. Conclusions

We have presented a framework that incorporates aleatoric uncertainty in the CNN. Doing so allows the model to learn another label for a given patch, which can be considered an additional degree of freedom in the solution for that patch. As expected, and confirmed in the results, the impact on the solution was greatest where the uncertainty values were large. These regions are typically found in areas that show strong changes in texture and tone, usually near the ice edge, and in regions marked by abrupt changes of lake ice concentration. Adding aleatoric loss and bathymetry information improved the accuracy of mapping water and ice (specifically black ice). Notwithstanding, model improvements in water mapping were highly visible in deeper zones whereas for ice detection, only the metrics for Zone 2 in Lake Erie improved.

A possible area to research is the extension of this present work to many other lakes in Canada monitored on a weekly basis by CIS and for which lake–wide ice fractions are determined. Our current methodology can be leveraged for operational use to provide ice maps with inherent uncertainties. It expands product usability and will make researchers aware of aleatoric uncertainty when incorporating ice fractions in their physical/numerical retrieval models in the form of direct integration of observation error variance or as a quality control flag. Therefore, our proposed product enhances their retrieval accuracy and final–product uncertainty reports.