M A Veganzones

M Dalla Mura

M Dumont

I Zin

J Chanussot

IMPROVED SUBPIXEL MONITORING OF SEASONAL SNOW COVER: A CASE STUDY IN THE ALPS

Keywords: Snow coverage area, SCA, spectral unmixing, MODIS, Alps

The snow coverage area (SCA) is one of the most important parameters for cryospheric studies. The use of remote sensing imagery can complement field measurements by providing means to derive SCA with a high temporal frequency and covering large areas. Images acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS) are perhaps the most widely used data to retrieve SCA maps. Some MODIS derived algorithms are available for subpixel SCA estimation, as MODSCAG and MODImLab. Both algorithms make use of spectral unmixing techniques using a fixed set of snow, rocks and other materials spectra (endmembers). We aim to improve the performance of a modified version of MODIm-Lab algorithm by exploring advanced spectral unmixing techniques. Furthermore, we make use of endmember induction algorithms to obtain the endmembers from the data itself instead of using a fixed spectral library. We validate the proposed approach on a case study in the mountainous region of the Alps.

INTRODUCTION

The snow coverage area (SCA) is one of the most important parameters for cryospheric studies [START_REF] Salomonson | Estimating fractional snow cover from MODIS using the normalized difference snow index[END_REF]. The amount of snow fallen on the surface is of fundamental importance for studies in many fields such as hydrology (for the retrieval of the snow water equivalent, or flood forecasting), energy (estimations of the expected hydroelectric power), glaciology (for estimation of mass balances of glaciers) and climate change studies (detection of anomalies and trends in the evolution of the snow coverage over the years due to climate fluctuations). The estimation of the SCA is usually conducted using field measurements, either manual surveys or automatic weather stations on the ground. Although this approach leads to accurate measurements, these are limited to local scale and do not provide data with a large spatial coverage. Remote sensing imagery can complement these measurements by providing means to derive SCA with a high temporal frequency and covering large areas through the analysis of the acquired imagery.

Among all the remote sensing products, images acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS) are perhaps the most widely used data to retrieve SCA maps. MODIS is an imaging sensor aboard the Terra and Aqua satellites, which are viewing the entire Earths surface every one to two days, acquiring data in 36 spectral bands 1 . The MODIS snow products include the MODIS fractional snow cover based on spectral unmixing (MOD-SCAG) [START_REF] Thomas | Retrieval of subpixel snow covered area, grain size, and albedo from MODIS[END_REF]. The MODSCAG product assumes that the measured reflectance of each pixel is given by a linear mixture of a set of pure spectra (e.g., snow, ice, vegetation, and rocks), which are called endmembers. Although the MODSCAG product has proved to be more accurate than the methods based on the normalized difference snow index [START_REF] Hall | Accuracy of the MODIS snow products[END_REF][START_REF] Painter | Assessment of methods for mapping snow cover from MODIS[END_REF], its wide application is limited by the selection of the endmembers. Currently, these are fixed and do not comprehensively account for the different states of the snow surface (e.g., dirty snow is not accounted, [START_REF] Mary | Intercomparison of retrieval algorithms for the specific surface area of snow from near-infrared satellite data in mountainous terrain, and comparison with the output of a semi-distributed snowpack model[END_REF]). The MODImLab algorithm introduced in [START_REF] Sirguey | Seven years of snow cover monitoring with MODIS to model catchment discharge in New Zealand[END_REF], corrects the reflectance values by considering topographic and atmospheric effects, and obtains the SCA with the algorithm described in [START_REF] Sirguey | Improving MODIS spatial resolution for snow mapping using wavelet fusion and ARSIS concept[END_REF] which performs a spectral unmixing on a set of images at 250m by fusing the reflectance bands at 500m with the two at 250m. Authors in [START_REF] Mary | Intercomparison of retrieval algorithms for the specific surface area of snow from near-infrared satellite data in mountainous terrain, and comparison with the output of a semi-distributed snowpack model[END_REF] extended the work in [START_REF] Sirguey | Seven years of snow cover monitoring with MODIS to model catchment discharge in New Zealand[END_REF] by comparing different retrieval techniques for snow grain size.

Here, we extend the approach in [START_REF] Mary | Intercomparison of retrieval algorithms for the specific surface area of snow from near-infrared satellite data in mountainous terrain, and comparison with the output of a semi-distributed snowpack model[END_REF] by exploring the application of more advanced spectral unmixing techniques for increasing the precision in retrieving the SCA maps. We address the spectral variability issue by using a (partially) constrained least squares unmixing (CLSU) algorithm. Furthermore, we make use of endmember induction algorithms to induce the endmembers from the data instead of using a fixed spectral library. By adapting the endmembers to the data we expect to improve the unmixing results while still being capable of identifying the different snow conditions. Experimental validation will be provided through a case study in the mountainous region of the Alps.

The remainder of the paper is as follows. In Sec. 2, an overview of spectral unmixing is given. In Sec. 3, advanced spectral unmixing techniques and the induction of endmembers from data are presented. In Sec. 4, we give experimental validation of the improved subpixel monitoring of snow cover by a case study in the region of the Alps. Finally, we provide some conclusions in Sec. 5.

SPECTRAL UNMIXING

The Spectral Unmixing (SU) is the process by which a given pixel is decomposed in the spectral signatures of the materials that contains and their fractional abundances. The most common approach relies in the assumption that a pixel is formed by a linear combination of the spectral signatures of pure materials present in the sample (endmembers), normally corresponding to macroscopic objects in scene, plus some additive noise. This is the so-called Linear Mixing Model (LMM). Let E = [e 1 , . . . , e m ] be the pure endmember signatures, where each e i ∈ R q is a q-dimensional vector. Then, the LMM states that a pixel signature, r ∈ R q , is defined by the expression:

r = m i=1 e i φ i + n, ( 1
)
where φ is the m-dimensional vector of fractional per-pixel abundances and n is an independent additive noise component. Some physical constraints can be enforced regarding the fractional abundances. The Abundance Non-negative Constraint (ANC), φ i ≥ 0, ∀i, ensures that there are not materials with negative contributions. The Abundance Sum-to-one Constraint (ASC), m i=1 φ i = 1, makes the total contributions of the materials sum up to one. In this model, the pixels lie inside a simplex whose vertexes are the endmembers.

Given the endmembers, the spectral unmixing aims to estimate their fractional abundances for each pixel independently. This can be modelled as an inversion problem:

φ = arg min φ r - m i=1 e i φ i 2 , (2)
subject to ANC and ASC abundance constraints. A solution to (2) is given by the full constrained (FC) least squares, which in the spectral unmixing context is known as the FCleast squares unmixing (FCLSU).

ADVANCED SPECTRAL UNMIXING AND ENDMEMBERS INDUCTION

Spectral variability

In real scenarios, the LMM and the fractional abundances estimated by the FCLSU are seriously affected by spectral variabilities due to variable illumination and environmental, atmospheric, and temporal conditions. By ignoring these variations, errors are introduced and propagated through hyperspectral image analysis [START_REF] Zare | Endmember variability in hyperspectral analysis: Addressing spectral variability during spectral unmixing[END_REF]. A simple way to address spectral variations is by dropping the ASC constraint, and solving (2) by (partially) constrained least squares unmixing (CLSU).

In this model, variabilities can be modelled as a scaling factor, λ > 0, that affects each endmember at a given pixel:

r = m i=1 e i α i φ i + n. (3)
In this extended LMM model, the hyperspectral pixels lie inside the positive hypercone defined by the endmembers. The CLSU implements the extended LMM (3) by solving the following inversion problem, r = m i=1 e i a i + n, where the estimated weighting factors, a i , are ANC but not ASC. The estimated weighting factors incorporate the information from the spectral abundances, θ i and the scaling factors, α i , that is, a i = α i φ i . In order to retrieve the fractional abundances, it is possible to assume that the scaling factor is the same for all the endmembers, α i = α, ∀i, and estimate it as α = m i=1 a i . Then, the fractional abundances are obtained by normalizing the vector of weighting factors by the estimated scaling factor, φ i = a i /α.

Endmembers induction algorithms

Most of the times, the spectral signatures of the materials are unknown, and the set of endmembers must be built by either selecting spectral signatures from a spectral library, or by inducing them from the image itself. Both can be performed manually or in an automatic way. In order to automatically induce the set of endmembers from the image, the use of some endmember induction algorithm (EIA) is required. The hyperspectral literature features plenty of such algorithms [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF].

Among them, algorithms exploiting the geometrical characteristics of the LMM provide a simple and powerful approach to endmember induction. These algorithms look for a simplex set that contains the data. Therefore, finding the endmembers is equivalent to identifying the vertices of the simplex. The vertex component analysis (VCA) algorithm [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF], presents a great trade-off between performance and computational complexity. VCA iteratively projects data onto a direction orthogonal to the subspace spanned by the endmembers already determined. The new endmember signature corresponds to the extreme of the projection. The algorithm iterates until all endmembers are exhausted.

A STUDY CASE IN THE ALPS

We explored a case study with MODIS images in the mountainous region of the Alps. We followed the processing chain detailed in [START_REF] Mary | Intercomparison of retrieval algorithms for the specific surface area of snow from near-infrared satellite data in mountainous terrain, and comparison with the output of a semi-distributed snowpack model[END_REF] and compared the use of FCLSU with respect to CLSU given a spectral library of 8 endmembers (see Fig. ??). We also compared the use of an endmember induction algorithm before applying CLSU. In the latter case, we used the VCA algorithm. The MODIS acquisitions comprises five snow seasons: 2005-06, 06-07, 09-10, 10-11 and 11-12. After removing those images where the scene was fully covered by clouds, 854 acquisition days remained. Fig. 2 shows the average Spectral Angle Mapper (SAM), average ERGAS and average Q index reconstruction errors. These are common quality measures in unmixing literature, whereas for the former two, zero is the best value; while for the latter, one is the best. The CLSU outperforms the FCLSU in terms of reconstruction quality, and the use of CLSU with the set of endmembers induced by the VCA gives the best reconstructions. It can be noted that the errors in some acquisitions are particularly bad for the FCLSU and CLSU approaches. This is due to particularly bad climatological conditions (wide cloud covering). The EIA+CLSU approach overcomes this issue due to the adaptation of the induced endmembers to the particularities of the data. Fig. 3 shows the average fractional abundances obtained by FCLSU and CLSU approaches respectively. It is interesting to note that the FCLSU approach reports a higher abundance of glacier snow, while the CLSU approach reports old coarse and ice instead. The peaks of 100% abundance corresponds to the cloudy acquisition days.

CONCLUSIONS

We have preform a preliminary study of the possibilities and challenges of using a partially constrained unmixing instead of the conventional fully constrained one, in order to address issues related to spectral variabilities. We also considered the use of endmember induction algorithms instead of a spectral library. The CLSU and the EIA+CLSU approaches results in an unmixing that fits the data better than the conventional FCLSU approach. However, the fractional abundances reported by the CLSU considerably differ from the ones reported by the FCLSU approach, where the latter reports an abundance of glacier snow, while the former reports abundances of old coarse and ice. This aspect will be further investigated by considering snow cover maps obtained from the high spatial resolution SPOT sensor and on field ground truth measures obtained by MeteoFrance. Also, the meaning of the endmembers induced by the EIA will be considered.

Fig. 1 .

 1 Fig. 1. Spectral library of 8 endmembers used to unmix the data.

Fig. 2 .

 2 Fig. 2. Average reconstruction errors for the Arves dataset: (top) SAM, (middle) ERGAS, (bottom) Q index.

Fig. 3 .

 3 Fig. 3. Average reconstruction errors for the Arves dataset: (top) FCLSU, (bottom) CLSU.

http://modis.gsfc.nasa.gov

Acknowledges

This work has been supported by a grant from Labex OSUG@2020 (Investissements d'avenir ANR10 LABX56).


François Chaumette

Visual servoing

Keywords:

Related Concepts

-Visual navigation -Visual tracking -Robot vision -Pan/Tilt/Zoom camera control -Image-guided robotic surgery

Definition

.

Background

Visual servoing relies on techniques from image processing, computer vision, and control theory. It can be seen as sensor-based control from a vision sensor. An iteration of the control scheme consists of the following steps:

acquire an image; -extract some useful image measurements; -compute the current value of the visual features used as inputs of the control scheme; -compute the error between the current and the desired values of the visual features; -update the control outputs, which is usually the robot velocity, to regulate that error to zero, i.e. to minimize its norm.

For instance, for the first example depicted on Figure 1, the image processing part consists in extracting and tracking the center of gravity of the moving people, the visual features are composed of the two Cartesian coordinates of this center of gravity, and the control schemes computes the pan and tilt velocities so that the center of gravity is as near as possible of the image center despite the unknown motion of the people. In the second example where a camera mounted on a six degrees of freedom robot arm is considered, the image measurements are a set of segments that are tracked in the image sequence. From these measurements and the knowledge of the 3D object model, the pose from the camera to the object is estimated and used as visual features. The control scheme now computes the six components of the robot velocity so that this pose reaches a particular desired value corresponding to the object position depicted in blue on the images.

Fig. 1. Few images acquired during two visual servoing tasks: on the top, pedestrian tracking using a pan-tilt camera; on the bottom, controlling the 6 degrees of freedom of an eye-in-hand system so that an object appears at a particular position in the image (shown in blue).

Theory

Main if not all visual servoing tasks can be expressed as the regulation to zero of an error e(t) which is defined by e(t) = s(m(r(t)), a) -s * (t).

(

The parameters in (1) are defined as follows [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. The vector m(r(t)) is a set of image measurements (e.g. the image coordinates of interest points, or the area, the center of gravity and other geometric characteristics of an object). These image measurements depend on the pose r(t) between the camera and the environment. They are used to compute a vector s(m(r(t)), a) of visual features, in which a is a set of parameters that represent potential additional knowledge about the system (e.g. coarse camera intrinsic parameters or 3D model of objects). The vector s * (t) contains the desired value of the features, which can be either constant in the case of a fixed goal, or varying if the task consists in following a specified trajectory.

Visual servoing schemes mainly differ in the way that the visual features are designed. As represented on Figure 2, the two most classical approaches are named image-based visual servoing (IBVS), in which s consists of a set of 2D parameters that are directly expressed in the image [START_REF] Weiss | Dynamic sensor-based control of robots with visual feedback[END_REF][START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], and pose-based visual servoing (PBVS), in which s consists of a set of 3D parameters related to the pose between the camera and the target [START_REF] Weiss | Dynamic sensor-based control of robots with visual feedback[END_REF][START_REF] Wilson | Relative end-effector control using cartesian position-based visual servoing[END_REF]. In that case, the 3D parameters have to be estimated from the image measurements either through a pose estimation process using the knowledge of the 3D target model, or through a partial pose estimation process using the properties of the epipolar geometry between the current and the desired images, or finally through a triangulation process if a stereovision system is considered. Inside IBVS and PBVS approaches, many possibilities exist depending on the choice of the features. Each choice will induce a particular behavior of the system. There also exist hybrid approaches, named 2-1/2D visual servoing, which combine 2D and 3D parameters in s in order to benefit from the advantages of IBVS and PBVS while avoiding their respective drawbacks [START_REF] Malis | 2-1/2D visual servoing[END_REF].

s s R c R c s s R c R c
Fig. 2. If the goal it to move the camera from frame R c to the desired frame R c * , two main approaches are possible: IBVS on the left, where the features s and s * are expressed in the image, and PBVS on the right, where the features s and s * are related to the pose between the camera and the observed object.

The design of the control scheme is based on the link between the time variation of the features and the robot control inputs, which are usually the velocity of the robot joints q. This relation is given by ṡ

= J s q + ∂s ∂t (2)
where J s is the features Jacobian matrix, defined from the equation above similarly as the classical robot Jacobian. For an eye-in-hand system (see the left part of Figure 3), the term ∂s ∂t represents the time variation of s due to a potential object motion, while for an eye-to-hand system (see the right part of Figure 3) it represents the time variation of s due to a potential sensor motion.

As for the features Jacobian, in the eye-in-hand configuration, it can be decomposed as [START_REF] Chaumette | Visual servo control, Part II: Advanced approaches[END_REF] J

s = L s c V n J(q) (3)
Fig. 3. In visual servoing, the vision sensor can either be mounted on the robot (eye-in-hand configuration) or observing it (eye-to-hand configuration). For the same robot motion, the motion produced in the image will be opposite from one configuration to the other.

where -J(q) is the robot Jacobian such that v n = J(q) q where v n is the robot end effector velocity; -c V n is the spatial motion transform matrix from the vision sensor to the end effector. It is given by

c V n = c R n [ c t n ] × c R n 0 c R n (4)
where c R n and c t n are respectively the rotation matrix and the translation vector between the sensor frame and the end effector frame, and where [ c t n ] × is the skew symmetric matrix associated to c t n . Matrix c V n is constant when the vision sensor is rigidly attached to the end effector, which is usually the case. Thanks to the robustness of closed loop control schemes, a coarse approximation of c R n and c t n is sufficient in practice to get an estimation of c V n . If needed, an accurate estimation is possible through classical hand-eye calibration methods. -L s is the interaction matrix of s defined such that s = L s v where v is the relative velocity between the camera and the environment.

In the eye-to-hand configuration, the features Jacobian J s is composed of [2]

J s = -L s c V f f V n J(q) (5)
where f V n is the spatial motion transform matrix from the robot reference frame to the end effector frame. It is known from the robot kinematics model.

c V f is the spatial motion transform matrix from the camera frame to the reference frame. It is constant as long as the camera does not move. In that case, similarly as for the eye-in-hand configuration, a coarse approximation of c R f and c t f is usually sufficient to get an estimation of c V f .

A lot of works have concerned the modeling of the visual features and the determination of the analytical form of the interaction matrix. To give just an example, in the case of an image point with normalized Cartesian coordinates x = (x, y) and whose 3D corresponding point has depth Z, its interaction matrix is given by [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]

L x = -1/Z 0 x/Z xy -(1 + x 2 ) y 0 -1/Z y/Z 1 + y 2 -xy -x (6)
where the three first columns contain the elements related to the three components of the translational velocity, and where the three last columns contain the elements related to the three components of the rotational velocity. By just changing the parameters representing the same image point, that is by using the cylindrical coordinates defined by ρ = x 2 + y 2 and θ = Arctan(y/x), the interaction matrix of these parameters has a completely different form [START_REF] Chaumette | Visual servo control, Part II: Advanced approaches[END_REF]:

L (ρ,θ) = -cos θ/Z -sin θ/Z ρ/Z (1 + ρ 2 ) sin θ -(1 + ρ 2 ) cos θ 0 sin θ/(ρZ) -cos θ/(ρZ) 0 cos θ/ρ sin θ/ρ -1 (7)
which implies that using the Cartesian coordinates or the cylindrical coordinates as visual features will induce a different behavior, that is a different robot trajectory and thus a different trajectory of the image point. Currently, the analytical form of the interaction matrix is available for most classical geometrical primitives, such as segments, straigth lines, ellipses, moments related to planar objects of any shape, and also coordinates of 3D points and pose parameters. Methods also exist to estimate off line or on line a numerical value of the interaction matrix. Omnidirectional vision sensors, the coupling between a camera and structured light, and even 2D echographic probes have also been studied. A large variety of visual features is thus available for many vision sensors.

Once the modeling step has been performed, the design of the control scheme can be quite simple. The most classical control scheme has the following form [START_REF] Chaumette | Visual servo control, Part II: Advanced approaches[END_REF]:

q = -λ J s + (s -s * ) + J s + ∂s * ∂t -J s + ∂s ∂t ( 8
)
where λ is a positivive gain tuning the rate of convergence of the system and

J s +
is the Moore-Penrose pseudo inverse of an approximation or an estimation of the features Jacobian. The exact value of all its elements is indeed generally unknown since it depends of the intrinsic and extrinsic camera parameters, as well as of some 3D parameters such as the depth of the point in Equations ( 6) and [START_REF] Hashimoto | Visual servoing: real-time control of robot manipulators based on visual sensory feedback[END_REF].

The second term of the control scheme anticipates for the variation of s * in the case of a non constant desired value. The third term compensates as much as possible a possible target motion in the eye-in-hand case and a possible camera motion in the eye-to-hand case. They are both null in the case of a fixed desired value and a motionless target or camera. They try to remove the tracking error in the other cases.

Following the Lyapunov theory, the stability of the system can be studied [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. Generally, visual servoing schemes can be demonstrated to be locally asymptotically stable (that is the robot will converge if it starts from a local neighborhood of the desired pose) if the errors introduced in J s are not too strong. Some particular visual servoing schemes can be demonstrated to be globally asymptotically stable (that is the robot will converge whatever its initial pose) under similar conditions.

Finally, when the visual features do not constrain all the robot degrees of freedom, it is possible to combine the visual task with supplementary tasks such as joint limits avoidance or the visibility constraint (to be sure that the target considered will always remain in the camera field of view). In that case, the redundancy framework can be applied and the error to be regulated to zero has the following form:

e = J s + (s -s * ) + (I -J s + J s ) e 2 (9)
where (I -J s + J s ) is a projection operator on the null space of the visual task so that the supplementary task e 2 will be achieved at best under the constraint that the visual task is realized [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]. A similar control scheme to ( 8) is now given by q = -λ e -∂e ∂t [START_REF] Malis | 2-1/2D visual servoing[END_REF] This scheme has for instance been applied for the first example depicted on Figure 4 where the rotational motion of the mobile robot is controlled by vision while its translational motion is controlled by the odometry to move at a constant velocity.

Applications

Potential applications of visual servoing are numerous. It can be used as soon as a vision sensor is available and a task is assigned to a dynamic system to control its motion. A non exhaustive list of examples are (see Figure 4):

the control of a pan-tilt-zoom camera, as illustrated in Figure 1 for the pantilt case; -grasping using a robot arm; -locomotion and dextrous manipulation with a humanoid robot; -micro or nano manipulation of MEMS or biological cells; -pipe inspection by an underwater autonomous vehicle; -autonomous navigation of a mobile robot in indoor or outdoor environment; -aircraft landing; -autonomous satellite rendezvous; -biopsy using ultrasound probes or heart motion compensation in medical robotics. -virtual cinematography in animation.

Fig. 4 .

 4 Fig. 4. Few applications of visual servoing: navigation of a mobile robot to follow a wall using an omnidirectional vision sensor (top line), grasping a ball with a humanoid robot (middle line), assembly of MEMS and film of a dialogue within the constraints of a script in animation (bottom line).