2.1 The Multi-Angle Snowflake Camera

The MASC is a ground-based instrument which automatically takes high-resolution and stereoscopic photographs of hydrometeors in free fall while measuring their fall velocity. Its working mechanism is only summarized hereafter, as more details and explanations can be found in Garrett et al. (2012), who provide an extensive description of the instrument. Three high-resolution cameras (2448 pixels × 2048 pixels), separated by an angle of 36^∘, are attached to a ring structure and form altogether the imaging unit (see Fig. 1). The focal point is located inside the ring at about 10 cm from each camera (with a focal length of 12.5 mm). Particles falling through the ring and detected by the two horizontally aligned near-infrared emitter–receiver arrays trigger the three flashes and the three cameras. The cameras' apertures and exposure times were adjusted in order to maximize the contrast on hydrometeor photographs while preventing motion blur effects, leading to a resolution of about 33.5 µm and a sampling area of about 8.3 cm² (see Praz et al., 2017). The maximum frequency of triggering is 3 Hz, which is three image triplets per second (see Fig. 6).

Figure 1(a) Side view of the MASC with the three flash lamps in white on top and the two detectors as white boxes on the side of the metal ring (in black and red in front). (b) Top view of the inside of the MASC, with the three cameras clearly visible.

Download

These specifications can be compared to the snow particle counter (SPC) which has been used in many studies of blowing snow (e.g., Nishimura and Nemoto, 2005; Gordon and Taylor, 2009; Kinar and Pomeroy, 2015; Guyomarc'h et al., 2019) and can be considered as the reference instrument for monitoring blowing snow (e.g., Crivelli et al., 2016). The SPC has a control volume of $\mathrm{2}\mathrm{mm}\times \mathrm{25}\mathrm{mm}\times \mathrm{0.5}\mathrm{mm}$ and assigns particles into 32 diameter classes between 50 and 500 µm. It provides information on particle diameter (assuming a spherical shape), particle number and particle mass flux usually at a 1 s resolution (but raw data are measured at up to 150 kHz, Nishimura et al., 2014). For more information about the SPC, the reader is referred to the articles mentioned above.

2.2 Data sets

The MASC data used to implement and validate the present algorithm were collected during three field campaigns. The first one took place in Davos, Switzerland, from October 2015 to June 2016. The MASC was placed at 2540 m a.s.l in a double-fence intercomparison reference (DFIR; see Fig. 2a), designed to limit the adverse effect of wind on the measuring instruments in its center (Goodison et al., 1998). The MASC was about 3 m above ground. The two other campaigns took place at the French Antarctic Dumont d'Urville (DDU) station, on the coast of Adélie Land, from November 2015 to February 2015 and from January to July 2017, in the framework of the Antarctic Precipitation, Remote Sensing from Surface and Space project (http://apres3.osug.fr, last access: 28 January 2020) (Grazioli et al., 2017; Genthon et al., 2018). The instrument was deployed on a rooftop at about 3 m above ground (see Fig. 2b). A collocated weather station and a micro rain radar (MRR) were also installed. Nearly 3 million images were collected during these measurement campaigns altogether.

Figure 2Experimental setup conditions of the MASC in a DFIR near Davos (a) and on top of a container at Dumont d'Urville (b).

Download

From this large number of data, subsets of pure precipitation and pure blowing snow images were manually selected and further analyzed to choose relevant descriptors and fit a two-component GMM. The task of selecting a sufficient number of representative images for both classes turned out to be more complicated than expected, in particular for the Antarctic data set in which mixed images are very frequent. Gossart et al. (2017) used ceilometer data collected at the Neumayer (coastal) and Princess Elizabeth (inland) stations in East Antarctica to investigate blowing snow, and they suggest that more than 90 % of blowing snow occurs during synoptic events, usually combined with precipitation. For the sake of generalization, a large number of representative events were selected across the three campaigns. The goal was to cover a wide range of hydrometeors types and a wide range of snowfall intensities for the precipitation subset. Similarly, a wide range of wind speeds and concentrations were considered to build the blowing snow subset. From the campaigns in Antarctica, pure blowing snow and hydrometeor events were highlighted by comparing time series of MASC image frequency, wind speed and MRR-derived rain rate, as illustrated in Fig. 3. It was noticed that during strong blowing snow events, the number of images captured by the MASC was much larger than during precipitation events (more than one image per second; see Fig. 6). Potential pure blowing snow events were selected when the MASC image frequency and wind speed were higher than their respective median estimated over the whole campaign (to select relatively high values), and no precipitation was detected during the preceding hour. Only events for which these criteria applied for over an hour consecutively were kept. To highlight pure precipitation, the principle was the same but the criteria were an image frequency and a wind speed lower than the median as well as a MRR precipitation rate greater than zero. The MRR has a certain detection limit, so it was noticed that events selected as blowing snow could also occur during undetected light precipitation. As a result, images from all events were rapidly checked visually, and the campaign logbook was consulted to ensure that the selection was consistent and coherent. In both cases, some events had to be removed because of obvious mixing of blowing snow and hydrometeors.

Figure 3Time series and scatter plots of MASC image frequency, wind speed (measured at 10 m) and MRR-derived rain rate for the Antarctica 2015–2016 campaign. The gray shading indicates days during which time steps have been selected for the training set as blowing snow (dark gray) or precipitation (light gray). In the bottom scatter plots, the markers figure the selected blowing snow and precipitation time steps. Points on the x axis in the left scatter plot are potential candidates for pure blowing snow.

Download

As the MASC was deployed inside a DFIR in Davos, no blowing snow events were selected from this campaign. Although the DFIR is supposed to shelter the inner instruments from wind disturbances, we noticed that many images do not solely contain pure hydrometeors. From a webcam monitoring the instrumental setup, we noticed that the fresh snow accumulated on the edges and borders of the wooden structure of the DFIR was frequently blown away towards the sensor. To enlarge the precipitation subset, events with high snowfall rate but not affected by outliers of fresh wind-blown snow were added. Finally, some sparse images of obvious pure hydrometeors in the middle of mixed events were also included in the training set. In total, each subset contained 4263 images and, despite possible remaining (limited) uncertainty in the exact type of images, is assumed to be accurate and reliable enough to serve as reference for the evaluation of the proposed technique (see Fig. 8 and Sect. 4.2).

Table 1Campaigns and dates of selected events for the blowing snow (BS) and precipitation (P) subsets.

Download Print Version | Download XLSX

3.1 Particle detection

The MASC instrument and the collected images are described in Sect. 2.1. Although a single particle activates the cameras, many MASC pictures contain multiple particles distributed over the entire image, especially when blowing snow occurs. In fact, the number of particles appearing on a single image is a key characteristic to distinguish between precipitation and blowing snow. As a result, it was deemed essential to detect all particles in each image rather than the triggering one only (which is sometimes unidentifiable). A key challenge of this approach was to get rid of the noisy background. For this purpose, a median filter was used. The brightness of the background strongly depends on the luminosity at the instant of the picture, which varies according to the time of day and can change abruptly in partly cloudy conditions when the sun suddenly appears from behind a cloud. As a result, the median filter shows better performance to remove the background when systematically recomputed over a small number of consecutive images. Assuming that snow particles rarely appear at the exact same position in several consecutive images, the median filter was chosen to be computed over blocks of five images per camera angle. To ensure complete removal of the background when its brightness is greater that the corresponding median, a factor of 1.1 was applied to the filter. Finally, as some limited residual noise can still remain in the filtered image, a small detection threshold of 0.02 grayscale intensity was applied to isolate the snow particles. Masks of the sky and reflecting parts of the background (i.e., metallic plates) were created for each camera. The multiplication factor and detection threshold are increased in the regions delineated by the masks if the normal filtering leads locally to more pixels detected that one can expect from real particles. These steps are illustrated in Figs. 4 and 5. Issues in the filtering may occur if consecutive images are separated by a period of time during which the ambient luminosity has changed significantly (e.g., before/after the sunrise or sunset). An example is shown in Fig. B1 in Appendix B.

Figure 4Raw image, median filter, filtered image and final binary image for an example of blowing snow particles. The image size is 2448 pixels × 2048 pixels, corresponding to 82 mm×68.6 mm. Original MASC images are in gray shades, but the color scheme used here aims to enhance contrast and details for visual purposes.

Download

Figure 5Raw image, median filter, filtered image and final binary image for an example of a hydrometeor. The color scheme is used to enhance details for visual purposes.

Download

3.2 Feature extraction

Machine-learning algorithms require a set of variables, commonly called features or descriptors, upon which the classification is performed. Because of the fragmentation of ice crystals when hitting the snow surface (e.g., Schmidt, 1980; Comola et al., 2017), blowing snow is expected to be characterized by much smaller particle size and much higher particle concentration than snowfall (e.g., Budd, 1966; Budd et al., 1966; Nishimura and Nemoto, 2005; Naaim-Bouvet et al., 2014). In this study, various quantitative descriptors were therefore calculated according to four different categories: the number of particles and their spread across the image, the size of the particles, the geometry of the particles, and the frequency at which the images are taken.

Since it is difficult to exactly guess which descriptors are the most adequate to differentiate between blowing snow and precipitation images, an extensive collection of features was extracted from the blowing snow and precipitation subsets and compared. The selection of the most relevant ones is explained in the next section. As the classification is performed at the image level, we need features at the same level, and the information on the geometry and size of each detected particle in the considered image must hence be transformed into a single descriptor for that image. Consequently, quantiles ranging from 0 to 1 and moments from 1 to 10 were computed out of the distribution of the considered feature within the image. The image frequency is a descriptor independent of the content of the image and thus from the detection of particles. It is therefore not affected by potential image-processing issues. As each image comes with its attributed timestamp, the average number of images per minute was calculated with a moving window. The full list of all computed descriptors is displayed in Appendix A. The extraction of features was conducted with the MATLAB Image Processing Toolbox, in particular the function regionprops (https://ch.mathworks.com/help/images/ref/regionprops.html, last access: 28 January 2020).

4.1 Feature selection and transformation

Selecting a relevant set of features and avoiding redundancy is essential for accurate classification, regardless of the classification algorithm. For each of the four categories of descriptors previously mentioned, the most relevant one (according to the criterion explained below) was kept. The descriptor maximizing the interclusters-over-intraclusters distance described in Eq. (1) was selected. This quantity represents the distance between the mean of the blowing snow and precipitation distributions (μ_BS and μ_P respectively), normalized by the sum of their respective standard deviations (σ_BS and σ_P respectively).

For the features describing the number of detected particles and their spread across the image, the cumulative distance transform was kept. It represents the sum over each entry of the distance transform matrix (https://ch.mathworks.com/help/images/ref/bwdist.html, last access: 28 January 2020) of the binary image. The distance transform matrix has the same dimensions as the binary image and computes, for each pixel, the Euclidean distance to the nearest 1 element (i.e., the nearest particle). As a result, an image with many particles well distributed over its entire surface will have a low cumulative distance transform, while a single particle, even particularly large, will have a high value. This descriptor is more robust to image-processing issues than the raw number of particles, as illustrated in Fig. B2 in Appendix B.

Concerning the size distribution of the particles detected in an image, the quantile 0.7 of the maximum diameter was selected (because it has the highest S value among the different quantiles tested). The maximum diameter (D_max) represents the longest segment between two edges of a particle (see Praz et al., 2017, for more details). A logarithmic transformation of this feature was performed to make the distributions of the two classes more Gaussian. The minimum (i.e., quantile 0) squared fractal index showed the greatest S value (hence discrimination potential) among the features related to the particle geometry indices. The fractal index (FRAC) is defined according to the formula proposed by McGarigal and Marks (1995) in the context of landscape-pattern analysis. It was also more recently used to quantify stand structural complexity from terrestrial laser scans of forests (Ehbrecht et al., 2017).

Due to its different nature, the image frequency descriptor was selected by default, but it is worth noting that it has the highest S value (Eq. 1) among all descriptors (Table 2). The marginal distributions of the selected descriptors for the training set are shown in Fig. 6 to provide an idea of their respective magnitude and variability, as well as to illustrate their discrimination potential. As noted above, the image frequency is the most informative descriptor to distinguish blowing snow and precipitation.

In summary, four descriptor categories (related to particle size, particle geometry and particle distribution within the image, and image frequency) have been defined to distinguish images collected during blowing snow or snowfall, based on the expected differences in particle size and concentration between the two. A number of descriptors were estimated from each image by computing various quantiles and moments of the distributions of geometric properties of the particles in the considered image. One descriptor from each of the four categories defined above (listed in Table 2) was then selected to be further used for classification as the one maximizing the interclusters-over-intraclusters distance defined in Eq. (1).

Figure 6Histograms of selected descriptors for the training blowing snow and precipitation subsets.

Download

Table 2Selected features and corresponding S values.

Download Print Version | Download XLSX

4.2 Model fitting

The choice for the binary classification task was made on a Gaussian mixture model, an unsupervised learning technique that fits a mixture of multivariate Gaussian distributions to the data (see Murphy, 2012; McLachlan and Basford, 1988; Moerland, 2000, for more details). The mathematical description of a multivariate normal distribution is provided in Eq. (2).

where x is a Gaussian multivariate random variable of dimension D, μ its mean and Σ its covariance matrix, with ^T the transpose operator.

The choice of an unsupervised approach is based on several reasons. First, unsupervised methods do not depend upon labels. Hence, it is not required to ensure correct labeling of each image in the training set. As mentioned earlier, many images are composed of a mix of blowing snow and precipitation, and it is thus difficult to guarantee the objectivity of all given labels. Second, a clear separation observed between the two subsets would be statistically highly significant as no prior information is provided to the learning algorithm about the classes. Third, for low-dimension problems, unsupervised methods are sometimes less prone to overfitting and have a better potential for generalization. A main advantage of the GMM compared to other unsupervised methods is to provide posterior probabilities on the cluster assignments and thus allow for soft clustering (i.e., probabilistic assignment). In the context of the present study, this is absolutely relevant as there exists a whole continuum of in-between cases of mixed images. It should be noted that the descriptors were selected using a reference set (see previous section), but the clustering conducted by means of the GMM is itself unsupervised.

A two-component GMM with unshared full covariance matrices was thus fitted to the four-dimensional ($\mathit{x}=\mathit{\{}{f}_{\mathrm{1}},f_{\mathrm{2}},f_{\mathrm{3}},f_{\mathrm{4}}\mathit{\}}$, where f_i represents the four features listed in Table 2) data composed of the blowing snow and precipitation subsets. The MATLAB Statistics and Machine Learning Toolbox was used for this purpose, and the model parameters were estimated by maximum likelihood via the expectation–maximization (EM) algorithm (https://ch.mathworks.com/help/stats/gaussian-mixture-models-2.html, last access: 28 January 2020). The features were standardized before fitting the model. The mixing weights (or component proportions) were artificially set to 0.5 by randomly removing 80 data points from the training set and fitting again the GMM to have perfectly balanced classes. This step is essential as the model will then be used to classify new images (possibly from other campaigns). There are no reasons to give more weight to one component, as the relative proportion of blowing snow and precipitation images strongly depends on the campaign location. The posterior probabilities are computed using the Bayes rule (Murphy, 2012):

where z_i is a discrete latent variable taking the values $\mathrm{1},\mathrm{\dots },K$ and labeling the K Gaussian components. $P(z_i=k|{\mathit{x}}_i,\mathit{\theta })$ is the posterior probability that point i belongs to cluster k (also known as the “responsibility” of cluster k for point i). $P\left({\mathit{x}}_i\right|z_i=k,\mathit{\theta })$ corresponds to the density of component k at point i (i.e., $\mathcal{N}\left({\mathit{x}}_i\right|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k)$), and $P(z_i=k|\mathit{\theta })$ represents the mixing weight (also denoted π_k). Note that the π_k's are positive and sum to 1. θ refers to the fitted parameters of the mixture model {${\mathit{\mu }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\mu }}_k$, ${\mathbf{\Sigma }}_{\mathrm{1}},\mathrm{\dots },{\mathbf{\Sigma }}_K$, ${\mathit{\pi }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\pi }}_K$}. P(x_i|θ) is the marginal probability at point i, which is simply the weighted sum of all component densities:

As the concern of this study is on two components only, a more compact notation will be used for the rest of the article. The latent variable z will be replaced by k_P and k_BS to refer to the precipitation and blowing snow clusters, respectively. The term θ, which denotes the model parameters, will be left implicit. Assuming we are at first interested in performing some hard clustering (i.e., single label to a given image), an image will be classified as blowing snow if $P\left(k_{\mathrm{BS}}\right|{\mathit{x}}_i)>P(k_{\mathrm{P}}\left|{\mathit{x}}_i\right)$. That is to say, if the posterior probability of belonging to the blowing snow cluster is greater than 0.5, an image will be classified as such (because the posterior probabilities sum to 1). The model performance was assessed by simply labeling the data points according to its initial subset. An overall accuracy of 99.4 % and a Cohen kappa score of 98.8 % were achieved. The Cohen kappa statistic adjusts the accuracy by accounting for correct predictions occurring by chance (Byrt et al., 1993). These high values indicate a very good performance of the fitted GMM. Figure 7 presents the fitted Gaussian components as well as the reference values (not used in the fitting) for each of the six possible pairs of the four descriptors. It clearly illustrates the performance of the fitted GMM and the discriminative power of the descriptor related to image frequency.

To investigate the stability of the Gaussian components, the precipitation and blowing snow subsets were both randomly permuted and divided into 10 equal parts. Ten new training sets of a balanced amount of each subset were created, and a new GMM was fitted. Figure 8 shows on the top line the boxplots of the Gaussian component parameters μ_d and σ_d (i.e., diagonal entries of Σ) for each of the four dimensions. The boxplots show a limited variability for each feature (below 10 %), indicating a reasonable stability of the fitted parameters. In addition, the bottom line of Fig. 8 presents the learning curves, and their fast convergence to the same horizontal line when more than 30 % of the training set is used indicates a data set large enough for a reliable fitting of the GMM, without overfitting.

Figure 7GMM contours and data points projected on the 2-D planes. The colors correspond to the four entries of the confusion matrix. The predictions result from the clustering and the ground truth is the given labels.

Download

Figure 8(a, b) Stability of the parameters μ and σ (diagonal entries of Σ) for the two Gaussian components. The boxplots show the distributions of these parameters for each dimension, after fitting the GMM on a 10-fold random split of the training set. The feature number follows the order given in Table 2. (c) Learning curves for the fitted GMM, showing the evolution of the train and Cohen's kappa test as a function of the proportion of the training samples used. The shaded areas correspond to the 25–75th percentile range computed over 40 iterations of 70 %–30 % random train-test splitting, and the bold lines are the medians.

Download

4.3 Flag for mixed images

As mentioned earlier, an asset of using a GMM model is the posterior probabilistic information that could help estimate the degree of mixing of an image. Data points located close to the decision boundary in the multidimensional space are likely to be composed of a mix of blowing snow particles and hydrometeors. However, distributions of posterior probabilities computed over thousands of new images from entire campaigns appeared to be stretched out on both ends of the domain (i.e., close to 0 or 1), and not many images were present in between. This is probably due to the nature of the descriptors and the resulting shapes and relative positions of the Gaussian distributions. In order to investigate this issue, an additional set of images corresponding to mixed cases was built: it exhibited clear differences in the posterior probabilities with the pure blowing snow and pure precipitation subsets. This differentiation was however around 10⁻⁶ (or $\mathrm{1}-{\mathrm{10}}^{-\mathrm{6}}$), which is not so informative as such. Consequently, it was decided to define a new index similar to the posterior probability of belonging to the blowing snow component but more evenly distributed across the range ]0,1[. The new index uses the negative logarithm of the posterior probabilities multiplied by the marginal probability. Taking the log of Eq. (3) for k_BS, we have (the same applies for k_P)

Noting that the term P(x_i|k_BS) on the right-hand side is $\mathcal{N}\left({\mathit{x}}_i\right|{\mathit{\mu }}_{\mathrm{BS}},{\mathbf{\Sigma }}_{\mathrm{BS}})$, one can substitute Eq. (2) into the above expression, which yields

The quadratic term on the right-hand side is the Mahalanobis distance, which is a distance that uses a Σ⁻¹ norm. Hence, it represents the distance between point x_i and the center of the distribution, corrected for correlations and unequal variances in the feature space (De Maesschalck et al., 2000). The second term is related to the determinant of the covariance matrix and equals −3.94 for the blowing snow component and −2.59 for the precipitation one. The two last terms are constant and sum to 4.37 (the component proportions were set to 0.5 and D=4). The right side of Eq. (6) is also known as the quadratic discriminant function (QDF, Kimura et al., 1987), commonly noted g_k(x_i). The minus in front of the logarithm on the left side of Eq. (6) is used to return positive values and facilitate subsequent graphical interpretations. Note that the constant term $\frac{D}{\mathrm{2}}\log \left(\mathrm{2}\mathit{\pi }\right)$ is often removed, but in this case it ensures that g_k(x) is positive, even for a Mahalanobis distance of zero. Figure 9 displays a scatter plot of the quadratic discriminant values of both components for the whole training set. The proposed index is defined as the angle of the vector representing a data point on the scatter plot, normalized by $\frac{\mathit{\pi }}{\mathrm{2}}$. It is thus computed as follows:

This normalized angle is bounded in ]0,1[, with values close to 1 (respectively 0) indicating a strong membership of the considered image (and not the respective proportions within this image) to the blowing snow (respectively precipitation) clusters. It is closely related to the asymmetry of the Mahalanobis distances between a point x_i and the centers of the two Gaussian distributions but corrected by the term $\frac{\mathrm{1}}{\mathrm{2}}\log \left(\right|\mathrm{\Sigma }\left|\right)$, which is different for the two components. The advantage of using the index in this form, rather than deriving it from the Mahalanobis distances alone, is to respect the decision boundary given by the maximum a posteriori (MAP) rule. This means a posterior probability of 0.5 yields a ψ index of 0.5. Finally, quantiles 0.9 (ψ_P0.9) and 0.1 (ψ_BS0.1) of the ψ index distributions of the points classified as precipitation and blowing snow, respectively, are retained as thresholds to flag potential mixed images. The idea is to allow, for both classes, 10 % of the training set images to be flagged as mixed. This value is qualitatively supported by the distribution shown in Fig. 9. It can be changed by the user to be more (increasing it) or less (decreasing it) strict on the classification as pure blowing snow or pure precipitation, depending on the intended application.

To ease reading and interpretation, a mixing index λ_m is introduced by linearly rescaling between 0 and 1 the ψ index of the images flagged as mixed (i.e., λ_m is not defined for pure precipitation or pure blowing snow images):

The mixing index also respects the hard clustering assignment boundary at 0.5: λ_m>0.5 indicates that the image contains a mix of blowing snow and precipitation particles, but it is closer overall to blowing snow and vice versa. Images with a normalized angle outside the two mixed thresholds have an undefined (NaN) index of mixing and are considered as pure blowing snow particles or pure hydrometeors. Results are provided treating all images independently, but the ψ index can also be averaged among the three camera angles to provide a unique value per image identifier as well. The median of the range (max – min) covered by the ψ values from the three individual views is about 0.08 in Davos and 0.05 at Dumont d'Urville, indicating a limited variability between the three views.

In summary, the classification as a mixed case is based on the angle characterizing the considered MASC image in the 2-D space formed by the axis related to pure blowing snow on the one hand and the one related to pure precipitation on the other hand. A mixing index λ_m is finally computed by linearly rescaling the normalized angle over the range of values corresponding to mixed cases.

Figure 9(a) Scatter plot of the quadratic discriminant values of both components for the training set. (b) Probability density functions (PDF) of the normalized angle for the precipitation and blowing snow subsets and thresholds to identify mixed images.

Download