Hydrometeor classification from two-dimensional video disdrometer data

The first hydrometeor classification technique based on two-dimensional video disdrometer (2DVD) data is presented. The method provides an estimate of the dominant hydrometeor type falling over time intervals of 60 s during precipitation, using the statistical behavior of a set of particle descriptors as input, calculated for each particle image. The employed supervised algorithm is a support vector machine (SVM), trained over 60 s precipitation time steps labeled by visual inspection. In this way, eight dominant hydrometeor classes can be discriminated. The algorithm achieved high classification performances, with median overall accuracies (Cohen’sK) of 90 % (0.88), and with accuracies higher than 84 % for each hydrometeor class.


Introduction
The two-dimensional video disdrometer (Kruger and Krajewski, 2002), 2DVD hereafter, significantly improves the capability of ground observations to describe the microphysics and microstructure of precipitation both in the solid and the liquid phase.The system, based on simultaneous observations of falling objects with two orthogonally oriented cameras, has been used to characterize the relationships linking raindrop shape, size, and terminal velocity (e.g., Thurai and Bringi, 2005;Thurai et al., 2009).It has also been employed to validate weather radar rainfall estimates (Schuur et al., 2001;Thurai et al., 2008;Cao et al., 2008;Zhang et al., 2008).Regarding snowfall, the 2DVD has been used to derive the statistical properties of particle size distributions of winter storms (Brandes et al., 2007), to improve the radar retrieval of equivalent liquid precipitation (Huang et al., 2010), and to simulate radar observations from measured snowfall microstructure (Zhang et al., 2011).
In the present work we employ 2DVD measurements for the classification of hydrometeors, with a special focus on ice-phase precipitation.The expression "hydrometeor classification" refers to techniques that aim to retrieve qualitative information about the dominant hydrometeor type that characterizes the precipitation.Such information can then be used for risk assessment (hazardous hydrometeor identification, like hail), for parametrization and validation of numerical weather prediction (NWP) models (e.g., Xue et al., 2000), or as support for microphysical investigations (e.g., Houze, 1993;Schneebeli et al., 2013).Today, hydrometeor classification techniques are implemented for different types of measurements.Typical examples in remote sensing are algorithms designed for ground-based polarimetric weather radars (Straka et al., 2000;Dolan and Rutledge, 2009;Chandrasekar et al., 2013), or for airborne radars and lidars observing ice-phase clouds (e.g., Shupe, 2007;Delanoe et al., 2013).These sensors enable the sampling of large domains at high resolution on a short timescale, but their classification retrievals are indirect, constrained by numerical simulations, and difficult to validate extensively.In contrast, airborne particle probe imagers (e.g., Feind, 2008) allow for direct classification along aircraft flight paths but only (given the high cost of these platforms) during intensive measurement campaigns.Ground-based instruments sample precipitation directly on site (although on small sampling volumes), and could be used to classify hydrometeors, thus becoming a point reference for remote sensing retrievals.Only a few research works have been devoted to the implementation of classification schemes for such instruments, and their focus was mostly on mixed-phase precipitation (Yuter et al., 2006) or the exploration of the potential synergy between multiple sensors (Marzano et al., 2010).Some commercial disdrometers (e.g., PARSIVEL; Löffler- Mang and Joss, 2000), originally designed for rainfall studies, provide a basic estimation of the precipitation type associated with each measurement by making assumptions on fall velocity and equivalent rainfall intensity.
In this context the information provided by the 2DVD is of particular interest because a pair of two-dimensional views, together with fall velocity, is provided for each particle.Such features alone allow expert users to interpret the images and visually recognize in them specific hydrometeor types (e.g., Zhang et al., 2011).This suggests that automatic classification methods, based on training over visually interpreted (labeled) episodes, may be well suited to perform hydrometeor classification.Supervised classification algorithms, such as the support vector machine, (SVM; Boser et al., 1992), are used today to perform similar kinds of tasks.For example, such techniques have been used in land cover classification (Camps-Valls and Bruzzone, 2005), wind power forecasts (Foresti et al., 2011;Zeng and Qiao, 2011), and weather prediction (Sullivan, 2009).The SVM is a linear and binary supervised classifier that finds the optimal separations between observations belonging to different classes.These observations are defined by a set of numerical features, and the optimal separation is learned from a training set in which the association between input observation and output class is known.The SVM is able to handle high-dimensional inputs, is less prone to over-fitting issues than other supervised methods (Camps-Valls and Bruzzone, 2009), and has been shown to perform relatively well on the prediction of weather types (e.g., Elmore, 2010).Furthermore, the SVM allows for the retrieval of the most relevant input features driving the classification, and can rank them in order of importance, with the implementation of multiple kernel learning (SVM-MKL) techniques (Rakotomamonjy et al., 2008;Tuia et al., 2010).
In this paper we train an SVM model on 2DVD data in order to classify eight hydrometeor classes of the dominant type of precipitation during time intervals of length t.Aggregation over time intervals is conducted to reduce the computational cost, which may be excessive if each particle is individually considered.A relatively short t of 60 s is chosen to minimize the effect of mixing of separate hydrometeor types.Individual 2DVD images are summarized over t with a high-dimensional set of numerical features, constituting the necessary input for the SVM classifier.Data collected in the Swiss Alps; in the French Jura; and in the southern part of Ontario, Canada, are used to train and validate the model.The manuscript is structured as follows.Section 2 describes the experimental setup and the basic 2DVD data.Section 3 presents the hydrometeor classification model.Section 4 presents the main results and their quality assessment, while Sect. 5 provides examples of the outputs of the hydrometeor classification.Section 6 concludes the paper and lists some future perspectives.

Experiment locations
The 2DVD data employed in the experiments were collected during three distinct field campaigns, between September 2009 and March 2013.The first campaign took place from September 2009 until June 2011 in Davos, Switzerland: the 2DVD was deployed in the Swiss Alps, at an altitude of about 2500 m a.s.l.Data for a total of 1700 h of precipitation in liquid, mixed, and solid phase were collected during this time frame.The second campaign took place in Remoray, France, from December 2012 until March 2013, at an altitude of about 920 m, in the context of an experiment focused on melting hydrometeors.A total of 270 h of precipitation in solid, liquid, and mixed phase was collected in this experiment.The third complementary campaign includes about 200 h of data (mainly solid precipitation) collected by three 2DVD instruments between December 2011 and March 2012 in the framework of the Global Precipitation Measurement mission (GPM, http:// pmm.nasa.gov/precipitation-measurement-missions), in the Cold-season Precipitation Experiment (GCPEx) that took place in Ontario, Canada.

2DVD instrument and data pre-processing
The 2DVD working principle is extensively described in Kruger and Krajewski (2002).Here we briefly summarize the most relevant features of the instrument.Figure 1 illustrates the 2DVD measurement principle (see Fig. 3 of Kruger and Krajewski, 2002, for more details).
Two orthogonal light sources coupled with two (A and B) line-scanning cameras generate two stacked measurement planes of about 11 cm×11 cm.The planes are vertically separated by a distance of around 6.2-7 mm (the exact value is determined by mechanical calibration).The cameras capture the falling particles at a resolution of 512 pixels (0.2 mm) at 34 kHz, and the vertical distance between the measurement areas of cameras A and B enables the measurement of fall velocity.
The raw images need to be processed before being employed.This involves the filtering of unreasonable measurements, as well as the rematching of the measurements taken from camera A and B in order to ensure that both images actually refer to the same particle.Filtering and rematching of 2DVD images is based on the work of Hanesch (1999) and Huang et al. (2010).We followed their methods, but with a noteworthy modification.Those studies, which were interested in snowfall only, restricted the maximum fall velocity to 4 and 6 m s −1 , respectively.We increased this upper boundary to 14 m s −1 , large enough to include, with sufficient margins, the range of variation found in rain (e.g., Beard, 1976) and large graupel (List and Schemena, 1971).
Despite this filtering, some non-realistic particles can still be observed in the output.These particles appear as large objects, vertically oriented and elongated, as shown in Fig. 2. Because of these peculiarities, they are easily identified and excluded from the analysis presented in this paper.The exact nature of these artifacts is unknown, but their vertical orientation and dimension suggest that they may be associated with small-scale wind effects, melting, or dripping, causing some particles to reside for an anomalous amount of time in the measurement areas of the two cameras.The proportion of rejected particles is on average 3 %, and it ranges between 0.5 and 13 % per day.A few precipitation events required higher rejection rates.They were excluded from the analysis presented in this work.Two additional potential sources of uncertainty (whose magnitude is currently not known in snowfall) are the image distortion that can occur when the horizontal component of the falling velocity of the particles is significant, and the local winds generated by the geometry of the instrument.To date, image distortion can only be corrected in rain, and in particular only for raindrops that possess an axis of rotational symmetry (Schönhuber et al., 2008).In contrast, the winds induced by the instrument itself can lead to an underestimation of particles having lower density and dimension 1 .Further research, which is beyond the scope of this paper, is needed to develop correction schemes for snowfall measurements in order to compensate for these two potential issues.

From single particles to global features
Pairs of 2DVD A-B images are available for each particle falling in the measurement area.For the purpose of the present work, it is useful to summarize this large amount of information by choosing a set of relevant descriptors2 .Then, the statistical distributions of these descriptors in a time step t are used as input information for the hydrometeor classification.The descriptors chosen in this work are listed in Table 1, and can be divided into three groups.

Joint descriptors
Two descriptors are obtained by combining the views of cameras A and B. They are particle falling velocity v [m s −1 ], and This descriptor was originally developed for raindrops, for which volumes could be calculated accurately from the 2-D views.It can be extended to particles of any shape as a reference measure of particle size.In the present work, D e is calculated according to the formulation of Huang et al. (2010).

Particle size
Other descriptors can be computed separately for camera A and camera B. Figure 3 illustrates some of them.The apparent shaded areas A A,B and perimeters P A,B are readily available from the 2DVD measurements, while thicknesses T A,B and widths W A,B of each particle are defined with respect to a bounding box around the particle (Fig. 3).v, D e , A, P , T , and W together describe the particle bulk dimension and velocity.

Particle shape
Additional descriptors are computed to better characterize particle shape.They are dimensionless shape metrics commonly used in the analysis of land cover images for remote sensing (Jiao and Liu, 2012), adapted for use on 2DVD images: where A r A,B [mm 2 ] is the area of the bounding box calculated for image A (or B).PF A,B is the pixel fraction and compares the shaded area with the area of the bounding box.PF A,B is an index of compactness, as is the roundness index (ROUND A,B ), which compares the shaded area with a circular approximation.FORM A,B and square pixel metric SqP A,B are shape complexity indices based on the area-toperimeter ratio (they increase with decreasing complexity), while fractal dimension FD A,B and shape index SI a,b are indexes based on the perimeter-to-area ratio (they increase with increasing complexity).ELONG A,B evaluates the degree of elongation of the particles.
As introduced above, the feature vector used in the SVM model refers to the distribution of descriptors in a time step t.Let us consider a time step, t, during which N particles are recorded.We compute the mean, median, some quantiles (10, 25, 75, 90 %), and interquantiles (Q75-25, Q90-10) of each descriptor over the N particles available.Additionally, for descriptors 3 to 13 of Table 1, we compute the correlation coefficient between the measurements of camera A and B. This leads to a set of 203 features calculated per time step: 16 derived from camera combinations, 88 calculated separately for A and B (so 176 in total), and 11 correlation coefficients.We selected t to be 60 s, as a trade-off between representativeness and temporal resolution.Additionally, no statistics are computed if N is lower than 20 particles for a specific time step (Appendix A).The 88 features calculated separately for A and B have been verified to be consistent between each other, with biases generally lower than 10 %.This suggests that the information carried by a single camera is sufficient for these 88 features.Therefore we can define, for each valid time step, a final feature vector x containing 115 useful features by using only the 88 single features from one of the two cameras.

Hydrometeor classification
This section details the proposed supervised classification approach.First we define the hydrometeor classes, then we detail how a training set is obtained, and finally we present the classifier employed and its implementation to the available data set.

Hydrometeor classes and training set
The principle of supervised classification methods is to use a set of N train labeled observations (or a training set) to train a classifier that will learn how to interpret new unlabeled observations.In our case, we need to assign the appropriate dominant hydrometeor type to a selected population of time steps of length t.The 2DVD offers the possibility to visualize the actual hydrometeor images, and the supervision was therefore conducted manually, according to the judgement of trained operators.Two operators independently interpreted the images by visualizing particle shapes, velocities, and taking into account the on-site environmental conditions (time of the year, temperature).Additionally, for the data collected in Davos, X-band radar observations over the region were available (e.g., Schneebeli et al., 2013), thus providing contextual information about the structure of the precipitation and, in stratiform cases, about the altitude of the melting layer.
The visualization and pre-interpretation of a wide range of time steps led to the selection of eight hydrometeor classes to describe the possible precipitation types in the available data set.Figure 4 shows an example of a typical particle belonging to each class.The classes are small-particle-like (SP), dendrite-like (D), column-like (C), graupel-like (G), rimedparticle-like (RIM), aggregate-like (AG), melting-snow-like (MS), and rain (R).The "-like" is added to emphasize that this approach identifies the dominant type of hydrometeor, which does not necessarily imply that (i) there is only one type of hydrometeor in the considered time step and (ii) that all hydrometeors exhibit pristine shape and geometry.
The definitions for some hydrometeor classes require clarification.SP time steps refer to particles falling during icephase precipitation that, given their size and the resolution of the instrument (0.2 mm), do not allow for proper visual shape recognition.Small aggregates, as well as single ice crystals, can be assumed to belong to this class.RIM is observed when riming processes smooth the shapes of the hydrometeors and increase their fall speed, while G time steps refer to fully developed graupel in which the original shape of the rimed crystal is no longer recognizable.MS is observed when the instrument records precipitation within the melting layer, and in these time steps raindrops, snowflakes, and smoother snowflakes with larger fall speed coexist in a mixed phase.
The creation of the training set involved the inspection of all the particles within each time step, in order to retrieve the dominant particle type and to provide the appropriate label.Particular attention was paid to select time steps that were as pure as possible, for the subsequent training of the classifier.The training set employed in the present work includes N train = 400 time steps, each of them numerically characterized by the 115 components of the associated feature vector x defined in Sect.2.3.

Classification method
In this section, we present the classifier used, the SVM.Then, we briefly detail an extension of the SVM that allows for the retrieval of the importance of each input feature (or group of features) in the model: SVM-MKL.

SVM
The support vector machine (Boser et al., 1992;Scholkopf and Smola, 2001), also known as the large margin classifier, is a linear classifier which finds the best linear separation between samples belonging to two classes.In our case, samples are time steps i of length t, represented by a vector x i of d = 115 features, and the classes are the dominant types of hydrometeors, y i .The model is trained on known couples v=1 .The SVM finds the best linear separation, of type f (x) = w, x + b, for which all training samples are at least at a distance of 1 from the separating plane.In other words, for all training samples, f (x) must be greater or equal to 1.To differentiate between positive and negative examples, we also multiply this expression by 1 if the sample is of the positive class and by −1 if it is of the negative class (the two types of hydrometeors).To summarize, the constraint is y i ( w, x i + b) ≥ 1, ∀i ∈ N train .The strategy pursued by the SVM (more details in Boser et al., 1992) is to find the separation which maximizes the distance between the closest points of each class, which are also called support vectors.This distance is called the margin and is inversely proportional to the norm of the parameters vector, i.e., ||w|| 2 .In order to allow for some classification errors, we also introduce a term xi i , which is non-zero for samples classified wrongly.The margin maximization problem is the following one: Complexity of the function C is a parameter that controls the constraint of perfect classification: if we allow some errors (by keeping C low), the margin becomes larger, thus reducing the dependence of the final model on training samples, which may be noisy or issued from errors in the measurements.A too high value of C drastically increases the value of the cost function as soon as errors are made.In this case, the resulting model will achieve perfect classification of the training samples, but the risk of over-fitting the training data and achieving poor generalization in the validation phase is higher.
This optimization model is solved using Lagrangian multipliers α, which allow us to rewrite the problem as When the optimal solution of Eq. ( 9) is found (i.e., the vector of coefficients α), the label of an unknown sample x v is assigned on the basis of the sign of the decision function, i.e., its position with respect to the hyperplane (f (x) = 0): In the present formulation, it can be observed that the SVM is only a binary classifier.A number of strategies exist to reduce multiclass problems to binary problems, and in the present work the one-against-one rule was employed (Hastie and Tibshirani, 1998).The one-against-one rule builds as many binary classifiers as there are pairs of classes.Each classifier is therefore used to assign the time step to one of two possible classes.The time step is eventually classified into the class that received the most assignments.

Nonlinear SVM
The SVM, as has been discussed above, can only solve linear problems (it defines a linear hyperplane).However, there is an elegant solution to solve nonlinear problems.Let us go back to Eqs. ( 9) and ( 10): the solution of the optimization does not depend on the training samples themselves but instead only on the dot products between samples (see x i , x j in Eq. 9).In the same way, the prediction for a new sample only depends on its dot products with the training samples (see x i , x v in Eq. 10).Dot products are measures of similarity between the samples.To perform nonlinear classification, we need to find an estimate of their similarity in a projected space of higher dimension H, where linear separation becomes possible3 .To avoid explicitly defining the coordinates of the samples in the projected space, i.e., φ(x i ), we can use functions that, even if expressed with points in the original space, correspond to dot products in the projected space H: these functions are called kernels.Without entering mathematical details, which the interested reader can find in Scholkopf and Smola (2001), a kernel corresponds to a similarity function such that K(x i , x j ) = φ(x i ), φ(x j ) .This means that, for a given projection φ(•), the kernel computed from x i and x j will correspond to their similarity in the space H defined by φ(•).A classification which is linear in the projected space is nonlinear in the original space, as illustrated in Fig. 5.
In practice, in order to obtain a nonlinear classification with the SVM, we replace the dot products in Eqs. ( 9) and (10) by kernel functions K(x i , x j ) and K(x i , x v ), respectively.A classical kernel to obtain such a behavior is the radial basis function (RBF), which is computed as follows: The RBF kernel acts as a Gaussian similarity, which is maximal when considering the same samples (K(x i , x i ) = 1), and decreases jointly with the decrease of similarity between the samples.The bandwidth σ controls the steepness of the Gaussian bell.

SVM-MKL
Even if very successful, SVM remains a black-box model, in the sense that no information about the importance of the initial variables can be retrieved from its results.All operations are optimized in the projected space H: this means that, while it avoids computation of projection of the samples explicitly, it also prevents the assessment of the importance of the different variables involved.Recent research has offered a solution to this problem by introducing the concept of multiple kernel learning (MKL;Rakotomamonjy et al., 2008).
SVM-MKL builds on the so-called Mercer conditions, stating that a weighted sum of any positive definite function (a requirement for all kernel functions) is again definitely positive (Mercer, 1905).This means that we can design a valid kernel by a linear combination of M base kernels K m (x i , x j ), each one considering single time step features (in this case M = 115) or sets of time step features (in this case M < 115 and equals the number of groups of descriptors): d m is the weight attributed to each kernel K m and is a measure of the importance of this kernel in the combination, i.e., of the variables composing it.It usually sums up to 1.The SimpleMKL algorithm proposed in Rakotomamonjy et al. (2008) alternatively optimizes the weights and the SVM and at the same time retrieves the relative importance of each group (d m ), and the SVM model associated with the final weighted combination.
In our experiment, we used SimpleMKL to find the best combination of a series of RBF kernels K m , each one assigned to a set of features referring to the same particle descriptor (M = 13; see Table 1).As an example, K 1 takes into account the eight statistical features (Q10, Q25, Q50, Q75, Q90, IQ75-25, IQ90-10, mean) associated with the hydrometeor fall velocity v descriptor, K 2 the eight features associated with the equivolumetric diameter D e , and so on.2. The elements C(i, j ) contain the number of observations classified in the ith class, which in reality belong to the j th class.The diagonal contains the correct classifications.
Given the confusion matrix, the global performance of the classifier is quantified by the overall accuracy (OA), and Cohen's kappa (K): where S is the total number of classes, and N the total number of observations (in our case S = 8 and N = N * val ).K takes into account the correct prediction that might occur by chance, namely P est , and is a robust metric in the case of unbalanced classes.
Then, we look at the performances obtained within each class.For this purpose, we use where OA k is the accuracy of the kth class, and POD k and POFD k are respectively the associated probabilities of detection and false detection.

Evaluation of the quality of the training set
N train observations are available in total as a training set, and we must verify that this amount is sufficient for the present task.In other words, we want to assess here whether a larger N train would significantly improve the hydrometeor classification.To do so we proceeded as follows: (1) N train = 400 was initially randomly split into N * train = 300 and N * val = 100; (2) N * train was iteratively reduced in size, while the original N * val was kept for validation; (3) evaluation of the performance was conducted at each step; and (4) steps ( 1)-(3) were repeated with 200 realizations of the original split.
Figure 6 shows the evolution of K as a function of the number of training samples in the training set (N * train ).We can observe that N * train larger than 200 did not lead to significant improvements in terms of K, while when N * train was smaller than 100, the performances started to degrade sharply.These results suggest that the total available labeled samples (400) are sufficient for the present classification task.

Evaluation of the classification performances
For validation purposes, we now focus on 200 realizations of the case N * train = 300, N * val = 100.The classification achieved accurate global results, both in terms of OA and K.As shown in Table 3, K and OA mean values were 0.88 and 89 %, and in 90 % of the cases they took values higher  The classification performance associated with each hydrometeor class is summarized in Fig. 7.It can be observed that all the hydrometeor classes were identified with median OA k always greater than 84 %, POD k greater than 0.84, and POFD k lower than 0.16.Overall, rainfall (R) hydrometeor class achieved the best scores, together with columns (C).R hydrometeors showed a POD k equal to 1, meaning that errors for this class were uniquely false detection.In contrast, C hydrometeors showed low POFD k and OA k very close to 1, and the errors for this class were due mainly to missed detections, with POD k with median scores around 0.9.Graupel (G) was mostly affected by missed detections, and showed a relatively large interquantile spread for POD k , around the median value of 0.87.Small particles (SP) had the highest false detection rate, with median POFD k close to 0.15.Dendritic snow (D) exhibited the largest interquantile spreads, around otherwise satisfactory median values of 87 % (OA k ), 0.87 (POD k ), and 0.13 (POFD k ), followed by rimed particles (RIM) that exhibited a similar behavior, achieving higher scores for all the metrics.Aggregates (AG) and melting snow (MS) were both correctly predicted, with lower interquantile spread, median OA k larger than 88 %, POD k larger than 0.88, and POFD k lower than 0.12.
A last consideration concerns the choice of SVM as a classifier.Other methods are used to solve similar tasks in various fields of the environmental sciences, for example linear discriminant analysis (LDA) or neural networks (NN) (e.g., Goosaert and Alam, 2009;Robert et al., 2013).Comparison with these two methods showed that the proposed SVM scheme outperforms LDA by more than 25% and NN by more than 15% in terms of K.

Ranking of descriptors
The SimpleMKL algorithm was applied to learn the most relevant descriptors in the classification process, as explained in Sect.3.2.3.Referring to Eq. ( 12), it was observed that 5 groups of features out of the 13 (1 per descriptor, each including the 8 or 9 statistical features extracted from its distribution in t = 60 s) accounted for about 70 % (Fig. 8) of the total weights and therefore are considered hereafter as the most important ones.They are, in decreasing order of importance, pixel fraction PF, velocity v, equivolume diameter D e , form index FORM, and thickness T , with associated weights d m of 0.193, 0.181, 0.13, 0.112, and 0.098, respectively.This does not imply that the remaining eight descriptors were negligible in the classification process, but that we expect to find a more immediate and intuitive physical meaning in these five top-ranked ones.

Application on unlabeled data
This section presents some examples of the classification output, on data not included in the training set of the algorithm, and collected during the measurement campaign in Davos.

17 March 2011
A snowfall event occurring on 17 March 2011 is presented in Fig. 9.The air temperature recorded in the very close vicinity (less than 50 m away) was constantly below freezing (≈ −5 • C) through the entire event duration, and different ice-phase hydrometeors were identified in the time window shown here.Initially (07:00-09:00 UTC) precipitation was dominated by small particles (SP), followed by a phase of instability (09:00-10:00 UTC) characterized by sharp variations in the identified hydrometeor classes.During the next relatively stable phase (10:00-12:30 UTC), graupel (G) and larger rimed particles (RIM) were identified.Panels (b), (c), and (d) of Fig. 9 illustrate the behavior in time of the three top-ranked particle descriptors, namely pixel fraction PF, equivolume diameter D e , and fall velocity v.The median PF was around 0.7 during the entire event, indicating relatively high particle compactness.The median D e was initially below 1 mm (SP phase), and it increased to values between 1 and 2 mm in the latter part of the event characterized mostly by G and RIM classes.v exhibited the same trends as D e , and it increased when rimed particles and graupel were dominant.

12 January 2011
A different situation is depicted in Fig. 10, relative to a snowfall event recorded on 12 January 2011.In this case, for the time window shown (19:00-24:00 UTC), precipitation was dominated by aggregates (AG) and dendritic-shaped snow (D, at the end of the event).By comparing the present case with the one shown in Fig. 9, we observe a wider range of variation of particle sizes, with D e ranging between 0.5 and 8 mm (AG).Particle compactness was lower, with median PF below 0.7 throughout the event, and slightly lower for D than for AG.This is due to the higher geometrical complexity of aggregates and dendrites relative to small particles and graupel.The velocity v did not exhibit peculiar trends, and it remained around values of 1 m s −1 in

5 August 2010
The precipitation event that occurred on 5 August 2010 (Fig. 11) illustrates the transition between liquid-phase and ice-phase precipitation well.In the first part of the event (05:00-07:45 UTC) air temperature was around 4 • C, and it dropped to 0 • C in the second part of the event (07:45-08:00 UTC).After 08:00 UTC the air temperature stabilized again around 0 • C.These trends in temperature are directly reflected in the dominant hydrometeor types classified: initially rain (R), then melting snow (MS), and finally aggregates (AG).The rain was characterized by small D e and 2 ≤ v ≤ 5 m s −1 (i.e., light rain), which is larger than the typical velocities of ice-phase hydrometeors, and very high compactness, with a median PF around 0.9.In the transition from R to MS and AG, a clear and relatively smooth trend was observed for the three descriptors shown: v decreased to median values around 1 m s −1 ; the spread of D e increased; and the median PF dropped to 0.6 in the AG phase at the same time, as the geometrical complexity of falling hydrometeors increased.
Generally, the transition between R, MS, and AG was well captured in the large available data set.Figure 12 shows the relative number of classifications for each of these three types of hydrometeors as a function of the temperature.Please note that temperature has not been used as an input in the proposed system (Table 1).R always occurred at positive temperatures, MS maximum occurrence was between 2 and 1 • C, and AG around 0 and −1 • C (in agreement with Hobbs et al., 1974).These results give us confidence in the ability of the proposed technique to provide a meaningful classification.

Summary and conclusions
In this paper we presented a hydrometeor classification method based on the interpretation of 2DVD data.The classification, conducted with the SVM technique, uses as input the statistical behavior of a set of particle descriptors over time steps with length t = 60 s.The SVM was trained with 400 examples labeled by expert users, and outputs the dominant hydrometeor type within t.Additionally, an estimation of the relative descriptive importance of the input features is provided, which is of particular interest when higherlevel information on the particle characteristics is required.
Discrimination is performed between eight hydrometeor classes: small-particle-like, dendrite-like, column-like, graupel-like, rimed-particle-like, aggregate-like, meltingsnow-like, and rain.Evaluation of the classification performance was conducted both in global and class-specific terms.The classifier achieved accurate results, with median OA and K of 90 % and 0.88, respectively.Each of the classes were identified with a median accuracy exceeding 84 %.
Additionally, once trained, the classifier is fast enough to be potentially implemented in real time.
Three classification examples together with the time evolution of the top-ranked particle descriptors were used to illustrate the typical classification products in pure snowfall events and in the transition between snowfall and rainfall.Global hydrometeor-type behavior as well as small-scale fluctuations could be observed.The proposed classification of hydrometeors from the 2DVD measurements provides additional information that can help to better understand the microphysical processes characterizing ice-phase precipitation events.This work has the potential to be a starting point for ground-based quantitative evaluation of products derived from polarimetric weather radars.It can also be adapted and implemented to receive inputs from other particle imaging systems (one or two dimensional), either ground-based or airborne, provided that human interpretation can be carried out for the particle in the training set and that geometrical descriptors can be computed from the particle images.
The main limitation is that the current implementation provides bulk information over a given time step of length t, which is large enough to be statistically significant, but cannot provide estimation of hydrometeor mixtures over t.Future work will focus on the development of a particle-byparticle classification -which is more challenging in terms of computational requirements -that can lead to explicit quantification of hydrometeor mixtures.

Figure 2 .
Figure 2. Example of A-B views of a non-realistic particle that needs to be filtered.

Figure 3 .
Figure 3. Examples of particles descriptors of 2DVD images (camera A and B).On camera A: width (W A [mm]) and thickness (T A [mm]) of the bounding box enclosing the particle.On camera B: particle apparent perimeter (P A [mm]) and shaded area (A A [mm 2 ]).

Figure 4 .
Figure 4. Examples of particle images (two camera views: A left, B right) belonging to time steps dominated by a particular hydrometeor class.

Figure 5 .
Figure 5. Illustration of the nonlinear SVM.(a) A nonlinearly separable data set in the input space X , involving two classes (squares and circles).(b) Projection on a 3-D space H by the kernel K(x i , x j ) = x i , x j 2 .(c) Linear classification in the projected space H (filled dots are support vectors).(d) Corresponding nonlinear decision function in the original space X .Adapted from Volpi et al. (2013).
The evaluation of the accuracy of classification is conducted via different metrics.The available N train training observations are divided into two parts (N * train and N * val ).N * train observations are used as a training set to optimize the SVM parameters C and σ and to train the SVM, while the remaining N * val observations are kept for validation.A comparison is made between the SVM classification output {y * i } an 8×8 confusion matrix C, as shown in Table

Figure 6 .
Figure 6.Evolution of K [-] as a function of the training set size.The solid red line indicates while the blue and brown areas represent Q75-25 and Q90-10, respectively.The size of the training set was varied with a step of 1 between 300 and 20.Statistics are based on 200 realizations.

Figure 7 .
Figure 7. Bar plots of (a) OA k [%], (b) POD k [-], (c) POFD k associated with the eight hydrometeor classes undergoing classification.Statistics were calculated over 200 realizations of the SVM validation.

Figure 8 .
Figure 8. Weights d m of the 13 K m kernels, associated with the 13 particle descriptors used in the present study.

Figure 9 .
Figure 9. Snowfall event recorded on 17 March 2011.Time series of (a) dominant hydrometeor type as classified with the SVM and local ambient temperature [ • C], as measured by a closely located (distance ≤ 50 m) weather station; (b) D e [mm]; (c) pixel fraction of camera A PF A ; and (d) fall velocity v [m s −1 ].In panels (b), (c), and (d), black dots connected by the red solid line indicate the median value, while the shaded areas depict Q90-Q10 and Q75-Q25, respectively.

Figure 12 .
Figure 12.Distributions of the occurrence of AG, MS, and R as a function of air temperature.The distributions are obtained by aggregation of all the 2DVD measurements collected during the field experiments in Davos 2009-2011 and Remoray 2012-2013, and temperature data are given by closely located weather stations.

Table 1 .
List of descriptors chosen to describe the particles recorded.Descriptors 1 and 2 come from the combination of camera A and B; 3 to 6 describe particle size, and 7 to 13 particle shape.

Table 2 .
Example of a confusion matrix obtained during validation of the SVM classification for a validation set N * val of 100 observations.Correct classifications are situated on the diagonal, and misclassifications are in the off-diagonal entries.

Table 3 .
Mean values and relevant quantiles of K [-] and OA [%] calculated over 200 iterations of the SVM validation procedure.