Andres E Bedoya-Velásquez

Marcos Herreras-Giralda

Roberto Román

Matthias Wiegner

Sidonie Lefebvre

Carlos Toledano

Thierry Huet

Romain Ceolato

Ceilometer inversion method using water-vapor correction from co-located microwave radiometer for aerosol retrievals

Keywords: Aerosol inversion, Vaisala ceilometer, microwave radiometry, ceilometer data pre-processing

Introduction

Atmospheric aerosols play a crucial role in atmospheric dynamics and the energy balance of the Earth. The main impact of the aerosols relatedinteractions are: (i) the aerosol-radiation interaction (ARI), affecting the radiative fluxes of the Earth by absorbing and scattering solar and thermal radiation, and (ii) aerosol-cloud interaction (ACI) which are mainly associated to the modification of cloud properties and precipitation caused by aerosols ( [START_REF] Boucher | Clouds and aerosols[END_REF]).During the last decades, different active and passive remote sensors in synergistic operation have become a powerful strategy for the better determination of the atmospheric aerosol properties (optical and microphysical).

Previous works have shown that synergy between active remote sensors as lidar systems (light detection and ranging) and passive remote sensors, e.g. sunphotometers or microwave radiometer (MWR), allows to obtain advanced and vertically resolved aerosol properties ( [START_REF] Chaikovsky | Lidar-radiometer inversion code (liric) for the retrieval of vertical aerosol properties from combined lidar/radiometer data: development and distribution in earlinet[END_REF][START_REF] Lopatin | Enhancement of aerosol characterization using synergy of lidar and sun-photometer coincident observations: the garrlic algorithm[END_REF][START_REF] Benavent-Oltra | Comparative assessment of grasp algorithm for a dust event over granada (spain) during charmex-adrimed 2013[END_REF][START_REF] Benavent-Oltra | Different strategies to 28 retrieve aerosol properties at night-time with the grasp algorithm[END_REF]) and to study phenomena like aerosol hygroscopic growth ( [START_REF] Bedoya-Velásquez | Hygroscopic growth study in the framework of earlinet during the slope i campaign: synergy of remote sensing and in situ instrumentation[END_REF]) and the aerosol vertical dynamics using as proxy the planetary boundary layer (PBL) height ( [START_REF] De Arruda Moreira | Analyzing the turbulent planetary boundary layer by remote sensing systems: the doppler wind lidar, aerosol elastic lidar and microwave radiometer[END_REF]). The main drawbacks of these synergies are the cost of having the instrumentation operating together and also that most of the instruments are semi-automatic, which means that qualified human operation is frequently needed.

Ceilometers are low power single-wavelength lidar-based instruments which operate automatically, unattended and continuously. These instruments are commonly used for cloud base height determination and PBL studies, but recently, ceilometers have become an useful alternative for aerosol studies such as typical robust lidar instruments. These systems have been widely spread along the world with more than 1000 ceilometers installed over Europe, Asia and America. Currently, the COST Action ES1303 TOPROF (TOwards operational ground based PROFiling with ceilometers, doppler lidars and microwave radiometers for improving weather forecasts) has dedicated part of their interests on working in a better characterization of the ceilometer products and related uncertainties, and also E-PROFILE, a program of EU-METNET (EUropean METeorological services NETwork), is focused on the harmonization of ceilometer measurements and data provision across Europe, meaning that the interest in quantitative aerosol products from these instruments is increasing. In the last decade, ceilometers started to be used for long term studies of phenomena less investigated with remote sensors like aerosol hygroscopic growth ( [START_REF] Bedoya-Velásquez | Long-term aerosol optical hygroscopicity study at the actris sirta observatory: synergy between ceilometer and in situ measurements[END_REF]), to improve the forecasting models for example to predict fog events ( [START_REF] Haeffelin | Radiation fog formation alerts using attenuated backscatter power from automatic lidars and ceilometers[END_REF]), to retrieve profiles of aerosol properties ( [START_REF] Wiegner | What is the benefit of ceilometers for aerosol remote sensing? an answer from earlinet[END_REF][START_REF] Jin | Ceilometer calibration for retrieval of aerosol optical properties[END_REF][START_REF] Cazorla | Near-real-time processing of a ceilometer network assisted with sunphotometer data: monitoring a dust outbreak over the iberian peninsula[END_REF][START_REF] Román | Retrieval of aerosol profiles combining sunphotometer and ceilometer measurements in grasp code[END_REF] and to characterize them ( [START_REF] Herreras | Evaluation of retrieved aerosol extinction profiles using as reference the aerosol optical depth differences between various heights[END_REF] and [START_REF] Titos | Retrieval of aerosol properties from ceilometer and photometer measurements: long-term evaluation with in situ data and statistical analysis at montsec (southern pyrenees)[END_REF] ).

CHM15k ceilometer model (Lufft manufacturer) is widely used for aerosol inversion, mainly because it operates with a similar configuration as the commercial lidar systems, therefore the quality of the signals have been deeply studied and their capabilities are well known [START_REF] Cazorla | Near-real-time processing of a ceilometer network assisted with sunphotometer data: monitoring a dust outbreak over the iberian peninsula[END_REF][START_REF] Román | Retrieval of aerosol profiles combining sunphotometer and ceilometer measurements in grasp code[END_REF]. Other ceilometers used for the same end are the CL31 and CL51 models (Vaisala Inc.), but as it has been shown in [START_REF] Kotthaus | Recommendations for processing atmospheric attenuated backscatter profiles from vaisala cl31 ceilometers[END_REF], depending of the firmware and other features, Vaisala ceilometers present some drawbacks such as non-expected signal shapes and high electronic noises. In [START_REF] Marcos | Analysis of four years of ceilometerderived aerosol backscatter profiles in a coastal site of the western mediterranean[END_REF] is presented a new type of correction that improves the signal shape, named dark signal removal, making a substitution of the dark current measurements. In addition, as the emission line of the Vaisala CL-51 ceilometer is centered around 910 nm, water vapor absorption plays a critical role affecting the quality of the signal. [START_REF] Wiegner | Correction of water vapor absorption for aerosol remote sensing with ceilometers[END_REF] describes a methodology to make a water vapor correction of the signal using modelled water vapor absorption cross section and radiosondes for retrieving aerosol properties. After considering all this pre-processing, it is possible to use the ceilometer signal for aerosol retrieval using traditional methods such as the Klett's algorithm ( [START_REF] Klett | Lidar inversion with variable backscatter/extinction ratios[END_REF][START_REF] Cazorla | Near-real-time processing of a ceilometer network assisted with sunphotometer data: monitoring a dust outbreak over the iberian peninsula[END_REF]).

The main objective of this work is to present a new methodology, based on a modified Klett algorithm [START_REF] Cazorla | Near-real-time processing of a ceilometer network assisted with sunphotometer data: monitoring a dust outbreak over the iberian peninsula[END_REF], to retrieve optical aerosol properties from Vaisala CL51 ceilometers. To this end, a data pre-processing is required, including the suppression of the dark current noise (DCN), height-dependant background (BG) noises, and water vapor correction in near-real time by using a co-located microwave radiometer (MWR). The methodology allows to determine the error propagation when modelled atmospheric profiles are used for aerosol inversion products instead of using co-located measurements.

The manuscript is organized as follows. The site, instrumentation and data used are presented in Section 2, while Section 3 explain the applied methodology in terms of the signal pre-processing. Section 4 describes the Klett method to retrieve the aerosol profiles, and Section 5 shows the uncertainty in the retrievals caused by the use of water vapor derived from modelled radiosoundings instead of MWR data. Finally, Section 6 summarizes the main conclusions.

Site, Instrumentation and data availability

2.1. ONERA site ONERA, the French Aerospace Lab, is a research institute located in Toulouse, Occitanie, in the southern part of France (N: 43 34'12", E: 128'24").

The mission of the Optictronics department (DOTA) is to conduct studies and research in Optronics. These studies are conducted primarily for the benefit of the fields of Aeronautics, Space and Defence, but also for other fields such as security, environment, astronomy and medical imaging. MELOPEE Lab is a light-scattering lidar laboratory dedicated to the development of active remote-sensing instrument for light-scattering investigations. For this work, a ceilometer and a ground-based microwaver radiometer located on roof-top of the building were used for this experiment. Toulouse is a region with a humid subtropical climate dominated by Autan wind, which is a south-easterly wind from the Mediterranean. Due to the Garone river that divides the city between east and west crossing it from south to north, Toulouse presents rather high relative humidity (around 80%) almost all over the year.

The seasonal behaviour drives to have hot summers and cold winters.

Vaisala CL51 ceilometer

A CL51 Vaisala ceilometer, located at the ONERA site, have been used in this work. This is an active remote sensor that operates continuously (24/7) emitting pulsed laser radiation towards the atmosphere centered at 910 ± 10 nm. The backscattered radiation by the atmosphere is collected by a telescope in coaxial configuration, reducing the overlap height.Unless the ceilometer is a new instrument with less than a year to be installed and the fact that manufacturer provided an overlap information that assures full-overlap after 50m, we have considered our products above 250 m agl because this work is not focused in the near-field measurements. The detection system is based on an APD detector. The backscattered signal from the atmosphere is measured with spatial and temporal resolutions up to 10 m and 15 s respectively. More technical information can be found in [START_REF] Kotthaus | Recommendations for processing atmospheric attenuated backscatter profiles from vaisala cl31 ceilometers[END_REF]. The spectral range of the emitted and received light by the instrument is affected by atmospheric water vapor absorption, which has a direct impact on the recorded attenuated backscatter profile, which is the main product of the instrument.

Considering a single scattering approximation and assuming that lidar constant (K) and overlap (O(R)) are well known, the ceilometer retrieved signal based on the elastic lidar equation can be written as follows:

RCS(R) K 0 O(R) = U (R) = β(R) exp -2 R 0 α(r)dr (1)
where RCS(R) is the range corrected signal (the recorded signal divided by the square of range); K 0 is a constant that involves system characteristics (optics and electronics); O(R) is the overlap function referring to the geometrical probability of fully signal collection as a function of height; β is the backscatter coefficient. The double pathway atmospheric transmittance is defined as T 2 (R) = exp -2 R 0 α(r)dr where α is the atmospheric extinction coefficient. In a general formulation, β and α takes into account the contribution of the particles and molecules in the atmosphere [START_REF] Klett | Stable analytical inversion solution for processing lidar returns[END_REF][START_REF] Klett | Lidar inversion with variable backscatter/extinction ratios[END_REF].

The U (R) term is the main product retrieved from ceilometers, but in practice this coefficient cannot be always determined because K 0 and O(R) parameters are unknown most of the time. Due to this elastic lidar formulation, the wavelength dependency is omitted from the equations. The analytical solution of the Eq.1 can be obtained by using [START_REF] Klett | Lidar inversion with variable backscatter/extinction ratios[END_REF]) assuming a constant extincton-to-backscatter ratio, also known as lidar ratio (LR). Other atmospheric parameters are involved in the solution of the lidar equation, but those will be discussed in further sections. The data used in this work were measured continuously from October to November 2019 with temporal resolution about 36 s/profile and vertical resolution of 10 m (from 10 m to 15400 km agl). The instrument records 2400 measures per day of RCS profiles, but it also provides cloud base height, and other metadata.

RPG-HATPRO MWR

A ground-based microwave radiometer (RPG-HATPRO, Radiometer physics GmbH) is co-located to the mentioned ceilometer. MWR is considered as a passive remote sensor that performs measures unattended of the brightness temperatures of oxygen and water vapor in the atmosphere. The oxygen is measured in the K-band (51-58 GHz) and the water vapor in the V-band from 22 to 31 GHz with a radiometric resolution between 0.3 and 0.4 rms errors at 1.0 s integration time.

A previously trained neural network algorithm ( [START_REF] Rose | A network suitable microwave radiometer for operational monitoring of the cloudy atmosphere[END_REF]) is used to retrieve temperature and relative humidity (RH) and absolute humidity (AH) profiles. The temperature and RH profiles are provided at 92 height bins with variable vertical resolution and covering the first 10 km of the atmosphere.

Temperature and RH profiles performance (accuracy and precision) has been studied in previous studies using radiosondes as references and finding that the temperature accuracy and precision is up to 1± K and close to 6± 8 % for RH under cloud-free conditions [START_REF] Bedoya-Velásquez | Seasonal analysis of the atmosphere during five years by using microwave radiometry over a mid-latitude site[END_REF]; this pointed out the potential of these systems to retrieve atmospheric variables. Regarding the data availability for this study, the MWR provides up to 600 temperature, relative humidity and absolute humidity profiles per day with temporal resolution up to 2 min/profile. Temperature profiles had a composite format that combined the high spatial resolution of the atmospheric boundary layer profiles product with the standard temperature profiles.

CIMEL sun/sky photometer

A sun/sky photometer CIMEL CE318-N (Cimel Electronique S.A.S.)( [START_REF] Holben | Aeronet-a federated instrument network and data archive for aerosol characterization[END_REF])

is operating since 2013 at the south-east part of Toulouse (43.57 N, 1.37 E, at 160 m asl), up to 8 km away from the ceilometer in straight line. This instrument provides automatic measurements of sun and sky radiation at several wavelengths (340, 380, 440, 500, 675, 870, 940, 1020 and 1640 nm).These measurements are processed by AERONET in order to derive the optical and microphysical aerosol properties integrated in column [START_REF] Holben | Aeronet-a federated instrument network and data archive for aerosol characterization[END_REF]. The main AERONET product is the spectral aerosol optical depth (AOD λ ). In addition, the channel of 940 nm is used for retrieving the total column water vapor (or precipitable water vapor). AERONET also uses AOD λ and almucantar sky radiance measurements to retrieve and provide additional aerosol properties including volume size distribution, complex refraction index, and single scattering albedo at various wavelengths ( [START_REF] Dubovik | A flexible inversion algorithm for retrieval of aerosol optical properties from sun and sky radiance measurements[END_REF][START_REF] Dubovik | Variability of absorption and optical properties of key aerosol types observed in worldwide locations[END_REF][START_REF] Dubovik | Application of spheroid models to account for aerosol particle nonsphericity in remote sensing of desert dust[END_REF]).In this work, the AERONET AOD data used is the AERONET level 1.5 (cloud-screened) from AERONET version 3 ([27]), with an uncertainty lower than ± 0.01 for the wavelengths larger than 440 nm and below 0.02 for shorter wavelengths.

Hysplit GDAS meteorological data base

As a result of the computer analysis and forecast calculation performed at the centers for Environmental prediction (NCEP), it is possible to use an operational system so called Global Assimilation Data System (GDAS) for running the Hybrid Single Particle Lagrangian Integrated Trajectory Model (HYSPLIT)( [START_REF] Stein | Noaa's hysplit atmospheric transport and dispersion modeling system[END_REF]). One of the HYSPLIT modules allows to retrieve modelled radiosondes with different spatial and temporal resolutions. In this work, we have retrieved HYSPLIT profiles of temperature, relative humidity and pressure for ONERA location with spatial resolution of 0.5 and 3 h of temporal resolution. The total database of radiosoundings used for this study is 488 profiles (October and November 2019) from 0 to 10 km.

Signal pre-processing

The methodology developed in this paper involves the following 2 stages (Fig. 1): (i) the signal pre-processing: contain dark current (DC), background (BG) and water vapor correction by using real co-located atmospheric measurements and (ii) refers to the semi-automatic Klett algorithm.

Dark current correction

First, a ceilometer data pre-processing is performed. For signal noise correction, two types of signal are taken into account. The first one is linked with the electrical noise of the detectors which is so called dark current noise values around zero until 9 km agl (see Fig. 2b), but above this height the signal is fluctuating strongly describing a S-shape. From 9 to 11 km agl positive curve can be observed, then the signal decrease from 11 to 14 km agl describing a negative curve for finally increase again from 14 km to the end of the profile. To focus the attention in the shape of the DCN, we plotted the smoothed the DCN in Fig. 2b, observing that the noise is quite oscillating around zero for daytime (due to the light contamination) than night-time, but the shape remains. Once the DCN is characterized, it must be directly subtracted from the range corrected ceilometer raw data.

According to [START_REF] Kotthaus | Recommendations for processing atmospheric attenuated backscatter profiles from vaisala cl31 ceilometers[END_REF] results, this behavior can be expected for Vaisala ceilometers, however this analysis let us evidence the impact that DCN suppression will have on the RCS raw signal. The signal-to-noise ratio (SNR) of this instrument decrease above 4 km agl, which was checked for the raw signal (not shown) and range corrected (see Fig. 2c). The RCS raw signal shown in Figure 2c is illustrating the RCS raw shape before any DCN subtraction.

The DCN signals shown in Fig. 2b and the corrected ones presented in Fig. 2d were smoothed by applying a filter (i.e. spatial-averaging) just for centering the attention on the main shape of the RCS raw without all of these electrical fluctuations shown in Fig. 2ac, but these smoothed profiles were never involved in the calculations. The first difference between black and gray profiles shown in Fig. 2d is the minimization of the S-shape fluctuation and secondly the reduction/suppression of some near-range electronic fluctuations.

After DCN correction, it can be noticed that the signal is noise cleaned up to 4 km agl, and the S-shape oscillation of the RCS raw has been reduced. This S-shape of the Vaisala ceilometer signal is related to the opto-electronic system noises, but this analysis is beyond of the scope in this work.

Background noise

The second noise evaluated during the data pre-processing is the environmental light contamination of the ceilometer signal,which plays a role as a bias to the signal and it is well-known as background (BG) noise. This calculation was performed after DCN correction to have a signal with a significant reduction in the electronic fluctuations (i.e. minimization of the last kilometers oscillation). The BG is a constant value commonly consider as the mean value in the last kilometers of the lidar signals (i.e. 1 kilometer or more), but only if the signal present a constant noise. In our case, the signals were not so well behaved, unless that, we have reduced the oscillations considerably by suppressing DCN. As we shown in Fig. 3, the BG presented a height-dependency from 7 to 15.4 km agl, therefore the following ranges of the profile were evaluated: BG1 (7-9 km agl), BG2 (9-11 km agl), BG3 (11-13 km agl) and BG4 (13-15.4km agl). One positive aspect observed on this height-dependency is that BG noise at each range presented a normal distribution (not shown here). Therefore, each range-dependent BG mean calculated for each profile is a good candidate to suppress the BG noise.

Finally, we use the BG value that improved the dynamic range of the profiles (i.e. allows us to have positive signal in a far range), meaning that we systematically analyse each profile, then we performed the DCN correction and finally chose the best BG as it shown in Fig. 3. Summarizing up, this procedure will enhance the probability of having good RCS profile to invert until 6 km agl, and it also improves the signal-to-noise ratio (SNR) in the medium/far range. The fact that ceilometers present a low SNR above 7 km agl suggest that only big events such as clouds and big aerosol plumes can be properly detected.

The signal retrieved from the Vaisala ceilometers presents a challenge to be pre-processed in order to retrieve aerosol optical properties. One of the main challenges is addressed in [START_REF] Kotthaus | Recommendations for processing atmospheric attenuated backscatter profiles from vaisala cl31 ceilometers[END_REF], where it was showed that most of Vaisala systems have a positive or negative signal distortion associated to electronic noise fluctuations. The correction of this fluctuation is tackled by suppressing cosmetic offsets ( [START_REF] Kotthaus | Recommendations for processing atmospheric attenuated backscatter profiles from vaisala cl31 ceilometers[END_REF]) or dark corrections ( [START_REF] Marcos | Analysis of four years of ceilometerderived aerosol backscatter profiles in a coastal site of the western mediterranean[END_REF]). In this work, the DCN measurement and also the search of the best BG value assure that signal noise level and shape is improved. The corrected signal is defined in the manuscript as follows

RCS * (R) = RCS raw -[DCN + BG] (2)
where RCS raw is the raw range corrected signal, DCN is a 30 min averaged profile and BG is the mean value that fulfilled the criteria already explained.

In Fig. 3 is presented the whole scheme of the signal pre-processing, showing that good BG selection is not a trivial process for ceilometers. In the example case showed in the top of the panel of Fig. 3ad, it is performed the pre-processing on 11 th October 2019 at 7 UTC and in the bottom of the panel (Fig. 3e-h) at 10 UTC. In black dots it is represented the signal spatially-averaged each 500m together with the standard deviation. This representation has two aims, in one hand to evidence the shape of the RSC raw and how it improves after the DCN suppression for different profiles. It can be noticed after this correction that the oscillation of the averaged points is substantially reduced in the far and near range. On the other hand, we can notice in Fig. 3b-f that there is a remaining oscillating noise above 7 km agl, therefore is imperative to evaluate the height-dependency of BG noise. The Fig. 3 shows that the ceilometer BG noise is changing from one measurement to another, therefore the application of a systematic BG correction must be considered to each profile separately and checking which BG value improves the dynamic range. From the example case at 7 UTC, we use the BG4, while at 10 UTC the BG2 was the one that improved the signal. With this procedure, we are increasing the amount of signal available to invert.

Water vapor correction

As mentioned before, one of the CL51 drawbacks is the impact of the water vapor absorption on the laser emission line Previous works have demonstrated that this effect can be minimized. The water vapor correction method used in this paper is based on the one proposed by [START_REF] Wiegner | Correction of water vapor absorption for aerosol remote sensing with ceilometers[END_REF]. To perform this correction, absolute humidity and temperature profiles are used as input. In this paper, those profiles were measured from co-located MWR each 2 min from 0 to 10 km agl. These profiles were interpolated to the ceilometer vertical resolution. Regarding the temporal coincidence of the measurements, the temporal resolution of the ceilometer was degraded to the MWR. For the water vapor correction it is necessary to rewrite the transmittance term presented in Sec.2.2, splitting the transmittance term into the contributions of the molecules, particles and water vapor (T 2 m , T 2 p , T 2 w ). In this work, the transmittance term for molecules is calculated by using Rayleigh theory fed with atmospheric measurements (temperature and RH) from the MWR. Transmittance from particles can be determined as a result of the crosschecking procedure during the Klett inversion explained below. T 2 w can be defined as:

T 2 w = exp -2 R 0 α(r) w dr (3)
Where α(r) w = σ(r) w N (r) w , with N (r) w the water vapor number concentration and σ(r) w is the absorption cross section at the emitted wavelength.

The N (r) w will be calculated from atmospheric measurements of absolute humidity profiles like N (r) w =7.25 x 10 22 AH R w , where R w =0.462 Jg -1 K -1 .

σ(r) w is calculated following the results presented on [START_REF] Wiegner | Correction of water vapor absorption for aerosol remote sensing with ceilometers[END_REF] for the absorption cross section simulated profiles. In this context, we consider a Gaussian shape of the ceilometer emission spectrum centered at 910 nm with ∆λ=3.5

nm ( [START_REF] Wiegner | Correction of water vapor absorption for aerosol remote sensing with ceilometers[END_REF]). As water vapor decrease with height, the Eq.3 is solved for the first 10 km of the atmosphere, but results are only shown until 3 km agl. The water vapor corrected profile is named hereafter as RCS w = RCS * (R)/T 2 w .

Figure 4 shows an example case where the water vapor correction is applied following the explained procedure. The atmospheric transmittance due to water vapor molecules is calculated using the profiles of temperature and AH from MWR, and from them calculate the n w . In the example case shown in Fig. 4b (red line), it is possible to evaluate the bias between RCS * and RCS w , meaning that ceilometer signal tends to be overestimated if water vapor correction is not applied (Fig. 4b, gray line). The larger differences between water vapor corrected and no-corrected profiles are within the first 1.5 km agl meanly because i) tropospheric water vapor molecules are more abundant at these altitudes, and ii) the MWR spatial resolution is better until 2 km agl. This increase of the overestimation is seen on the bias plot (Fig. 4c), reaching almost 4 (a.u). The bias profile is highly noisy above 1.8 km agl which is caused by the increasing of the RCS * noise already discussed.

Retrieval of aerosol profiles

Ceilometer inversions have been successfully proved in Vaisala ceilometers in previous works ( [START_REF] Marcos | Analysis of four years of ceilometerderived aerosol backscatter profiles in a coastal site of the western mediterranean[END_REF]; [START_REF] Wiegner | Correction of water vapor absorption for aerosol remote sensing with ceilometers[END_REF]). However one of the main challenges remain in the use of ancillary information from models for correcting the signal and then invert it. In this work, we propose an approach to tackle this problem by combining co-located atmospheric profiles measured next to the ceilometer and then applying a modified methodology of the semi-automatic lidar Klett inversion proposed in [START_REF] Cazorla | Near-real-time processing of a ceilometer network assisted with sunphotometer data: monitoring a dust outbreak over the iberian peninsula[END_REF]. The aerosol retrieval method has been divided in 3-step process as follows:

• Signal smoothing (Step 1): First the signal is noise-cleaned applying a 1h average to the analyzed profiles. Then, in order to remove big peaks that remained from electronic noise, a spatial filter (40 m moving average filter) is applied. The aim of doing this filtering is to preserve the main shape of the RCS w with noise peak reduction.

• Rayleigh fit (Step 2): The Rayleigh fit procedure has been done by normalizing the signals, β m and RCS w , in order to have comparable magnitudes. After this normalization, it is possible to check automatically the region where both signals have similar slopes by means of a linear fitting between 2 to 5 km agl. This altitude range is pre-set as an input to the algorithm in order to apply the slope's method, starting from 2 km agl; however, in turbulent days (where the Atmospheric Boundary Layer(ABL) is too high), the presence of aerosol may remain at those altitudes, for that reason, we consider a range until 5 km agl, assuring that the slopes will be checked for different atmospheric volumes (each 60m) until the fulfillment of the criteria. The criteria to consider that both slopes are close enough is based on the R-squared is the correlation coefficient (the Goodness-of-Fit of the linear regression), then when r 2 > 0.7 between normalized β m and RCS w , and the percentual relative error between the slopes is below 0.2, the Rayleigh zone is identified. Once that a layer containing the Rayleigh information has been found, just one reference point inside of it has to be selected. The iterative methodology to select this Rayleigh reference point is described below

• Backward Klett inversion (Step 3).
The aerosol retrieval starts from the following lidar equation

RCS w (R) = Cβ(R)T 2 m T 2 p ( 4
)
Then the backward Klett equation can be expressed as follows,

β p = A 1 (R) A 2 (R) -β m (R) (5)
where

A 1 (R) = RCS w (R) exp 2 r ref r (LR p -LR m )β m dr (6)
and

A 2 (R) = RCS w (R) β(r ref ) m + β(r ref ) p + 2 r ref r LR(r) p A 1 (r)dr (7)
where β m is determined from Rayleigh theory using measured atmospheric profiles as inputs and the well-known molecular lidar ratio The results are shown in Fig. 5 for the case of 11 th October 2019 at ON-ERA site (occitanie Toulouse). The molecular profile (black dotted line in Fig. 5a) is calculated using Rayleigh theory fed with temperature and RH profiles from MWR. The cases presented here had a normalization range from 1.5 to 5 km, assuring a region where Rayleigh and RCS w slopes are comparable. In Fig. 5a the Rayleigh zone is detected from 1.8 to 3.5 km agl, and particularly in this example case a relative slope error of 0.05 with r 2 = 0.9 was found. The aim in the determination of the Rayleigh zone is to find the right altitude of reference to perform the backward Klett inversion, and as we see in Fig. 5, due to the oscillations that remain in the RCS w , the algorithm may found a valid reference point in the far ranges within the zone selected (i.e. above 3 km agl as it is shown in the figure). For that reason α p profiles are showing some variability in the proximity to the reference point(up to 0.006km -1 ).

In order to assure the full overlap, it is used 250 m agl as the minimum height during the Klett inversion. The example case shown here presents the evolution of the inversions since early morning until afternoon. The dynamic of the atmosphere started with the aerosol relatively compressed below the ABL (below 0.8 km agl) due to the lack of solar radiation, but after 10 UTC the layers started to mix along the atmospheric column showing two aerosol accumulations, from 0.4 to 0.5 km agl and from 1.0 to 1.2 km agl. According to AERONET, along the day, we had the presence of mixed aerosols between urban and dust particles over Toulouse, predominating the coarse mode in the size distribution reaching 0.019 µ m 3 /µ m 2 . The iterative Klett procedure for having β p profiles was performed by following the methodology, and then α p profiles were obtained as it is shown in Fig. 5b. During the morning, α p was up to 0.03 km -1 in the first atmospheric layers (from 0.2 to 0.5 km agl), and then from 10 to 15 UTC α p reached up 0.05 km -1 below the first 0.5 km agl. α p =0.03 km -1 for the pronounced aerosol peaks from 12 to 15 UTC.

One of the pursued parameters on lidar inversion is the LR, in our case, over Toulouse area there is no previous information about this parameter, therefore as one of the main goals for the application of a semi-automatic Klett is to have in the near future a robust data base that will help us to understand better the aerosol behaviour over the region and also to have more tools for retrieving aerosol optical properties. Table1 presents the LR obtained from Klett inversion, reporting an increasing LR along the day ranging from 58.1 to 99.0 sr, with a positive correlation with the AOD increase. In the early morning, the LR values are comparable with those reported in [START_REF] Cazorla | Near-real-time processing of a ceilometer network assisted with sunphotometer data: monitoring a dust outbreak over the iberian peninsula[END_REF] and [START_REF] Marcos | Analysis of four years of ceilometerderived aerosol backscatter profiles in a coastal site of the western mediterranean[END_REF] for dust aerosols up to 60 Sr, but after 10 UTC the LR values obtained here increased up to 99 sr. In most of the studies, the LR is imposed for the klett calculation because of the previous knowledge of the aerosol type, but in our case LR is obtained from the continuous AOD cross-checking during the algorithm iterations, giving us the possibility to explore a wide range of LR values that are related with different aerosol types. The larger values found for LR with Klett inversion might be associated to the presence of the mixed aerosol (e.g. polluted dust), the physical separation between the sun photometer and the ceilometer, and also due to instrumental parameters that increase the uncertainties (e.g lidar constant).

The inversions left two open paths for future works, from a technical sight will be the calculation of the ceilometer calibration constant along the seasons to considering different scenarios, and the design of the field campaigns for co-located measurements between ground base lidars and the ceilometer in order to inter-compare inversion products and improve the overlap information of the instrument. From a scientific point of view, we are studying the sensibility of the algorithm applied on synthetic lidar data and vaisala ceilometer signals in order to study the error propagation linked with the Rayleigh reference and LR values obtained along the iterations, and their impact over aerosol properties retrieved.

Errors from water vapor uncertainties

The methodology applied to retrieve aerosol profiles from vaisala CL51 ceilometer needs the a signal deep noise pre-processing and water vapor correction as we have shown in this work.Here, we have used atmospheric variables from the co-located MWR to improve the results, however, the water vapor correction in many stations with CL51 ceilometer is performed by taking atmospheric information from models, thus an extra uncertainty is added to the RCS and then propagated to the products retrieved. For that reason, we dedicate this section to quantify the errors when modelled inputs are considered instead of measured ones. The first part of this section is devoted to quantify the errors between atmospheric profiles calculated from HYSPLIT model instead of those measured with the MWR. In order to statistically characterize the errors between the model and MWR, we use the mean root mean square error (RMSE), and the mean bias error (MBE). These statistical quantities were calculated at each altitude in order to evaluate how far are the modelled data from the measured ones. The second part of this section is focused on the error quantification and its impact over theRCS w , β p and α p under both scenarios, assuming modelled data and measured ones. For doing that, the percentual relative error has been used (RE), where MWR quantities are considered as the reference.

HYSPLIT vs MWR data

In order to characterize the inversion calculation, a 2-step error estimation procedure is proposed (see Fig. 1). In the first step, it was used a 2 months database (October-November 2019) from atmospheric measurements of temperature and AH from MWR, and modelled radiosondes retrieved using the GDAS meteorological database from HYSPLIT model. The aim is to estimate the root mean square error (RMSE) and mean bias error (MBE) between temperature, AH and N (r) w profiles assuming MWR as the reference, evaluating the error between model and MWR atmospheric data.

Figure 6 presents the result of the two-months error comparison between HYSPLIT modelled radiosondes retrieved at ONERA location using the GDAS meteorological database and the MWR profiles. The data were chosen assuring coincidence of the temporal and spatial resolutions between soundings and MWR profiles from 0 to 10 km agl. For the analysis, 488 coincident cases were found. Temperature and AH, the main atmospheric variables involved in the pre-processing and inversion, were evaluated. N w was also analyzed but the RMSE and MBE errors were relatively lower and highly dependant of the temperature and AH errors. but after that altitude the RMSE increase monotonically until 6 • C. Therefore assuming the MWR as the reference, the results indicate that modelled radiosondes can be far from the MWR measurements at least 4 • C in the crucial zone for ceilometer profiles. The Fig. 6b shows that temperature retrieved profiles from HYSPLIT are always sub-estimating the MWR ones, being critical around the first hundred meters agl. For AH profiles (blue lines, Fig. 6a,b), the RMSE shows values below 1.5 gm -3 for the whole profile, however between 3.7 to 5.7 km agl, the RMSE presented a peak, which might be associated to the atmospheric region where the number of water vapor molecules decrease with height drastically as it was also seen in [START_REF] Marcos | Analysis of four years of ceilometerderived aerosol backscatter profiles in a coastal site of the western mediterranean[END_REF].

The MBE of the AH is quite variable, the ranges that overestimate the profile are from 0 to 200 m and 3 km to 5.8 km agl and the sub estimation ranges are from 1 km to 2.3 km and from 5.8 km to the end.

Uncertainty propagation on aerosol profiles

The second part of the error calculation is performed just to six particular cases where lidar inversion was applied. In the last section, it was seen the performance of the atmospheric variables between measured and modelled data, but here we centered our attention in determining the error propagation in terms of the percentual relative error (RE) for temperature, AH, RCS w , β p and α p in order to evaluate the error propagation during inversion process considering modelled data instead of atmospheric measurements.

In Fig. 7 are presented the errors of the input measurements used for retrieving β p and α p applying the Klett inversion already discussed. In Fig. 7a is reported the percentual RE between temperature and AH profiles considering MWR as reference, from there it can be seen that RE in temperature (red line) are higher in the first 800 m agl (up to 5 %), but then, the RE do not surpasses 1 %. Meanwhile, AH presented higher RE. AH relative errors (blue line) below 2 km agl were lower than 20 %, but above 2 km agl the errors can reach 100 % faster. The absolute errors on N w at using the model or the MWR data are really lower, practically negligible (not shown here).

Figure7b shows the RE on the RCS w profile for the same example case presented in Sec.4. Therefore, the results of having up to 1 % RE in temperature and close to 20 % in AH are causing an increase in the RCS w RE up to 30 % below 2 km agl, and above this height the absolute error increase drastically reaching up more than 100 % above 3.5 km agl. The evaluation atmospheric profiles for improve the main vaisala product. Additionally, the RE found here for RCS w have to be added to the high noise of this product, which in fact leads to increase the error peaks shown in Fig. 7b Table2 reports the mean values of the RE in percentage for six cases evaluated during the two-months study. The cases 1,2,3 were retrieved on 11 th October 2019 at 06:00, 09:00, and 12:00 UTC, case 4 to 20 th November Table 2: Mean RE calculated for six inversion cases using the semi-automatic Klett method. The variables evaluated were Temperature, AH, RCS w , β p and α p . The study cases 1,2,3 were measured on 11 th October 2019 at 6 h, 9 h, and 12 h respectively, case 4 on 20 th November 2019 at 6 h and cases 5 and 6 on 21 st November 2019 at 9 h and 12 h. The mean values are calculated for five different atmospheric volumes: from 0 to 0.5 km agl, from 0.5 to 1 km agl, from 1 to 2 km agl 0 to 0.5 km agl 0.5 to 1 km agl 1 to 2 km agl T (%) Case 1,8 ± 0,9 0,1 ± 0,1 0,70 ± 0,08 Case 3,9 ± 0,7 1,1 ± 0,3 0,65 ± 0,02 Case 3,4 ± 0,4 2,5 ± 0,4 1,3 ± 0,2 Case 1,5 ± 0,5 0,14 ± 0,08 0,5 ± 0,3 Case 2,1 ± 0,7 0,3 ± 0,1 0,6 ± 0,2 Case 2,0 ± 0,4 increase of the noise that we already discussed in the pre-processing section, linked with the quality of the ceilometer signals and also in Fig. 7b where the RE of the RCS w increased considerably faster with height, pointing out that aerosol inversion with ceilometers above 3 km agl might lead to larger errors.

conclusions

The purpose of this work is to continue tackling the central problems that Vaisala CL51 ceilometers have for improving aerosol inversion products.

Firstly, the ceilometer pre-processing is discussed by taking advantage of the termination hood external tool. Between day and night the instrument presented only small differences in terms of shape and noise levels of DCN signals detected, therefore we decided to work with 30 min-averaged nighttime DCN profiles. For the BG analysis, we used a systematic selection of the best height-dependence BG, avoiding that RCS * goes to negative values within the first 5 km agl. After suppressing DCN and BG, the signal remained positive for more than 5 km agl and S-shape oscillations were minimized.

The methodology designed involved the water vapor correction by using a co-located MWR for measuring atmospheric variables. This synergy allowed us to improve the quality of the signal by calculating the water vapor transmission term of the lidar equation using measurements of temperature and AH, and deriving from them the N w . For the inversion, Rayleigh profiles were calculated using MWR measurements and the RCS w profiles.In addition, the sun photometer AOD time series were interpolated to 910 nm in order to compare the results with those obtained from the integrated α profiles retrieved after Klett's inversion procedure. For the iterative Klett inversion method, a minimum difference between AOD 910 and integrated α profiles lower than 0.001 was set as the condition to constrain the aerosol amount along the iterations.

An error propagation was performed applying the methodology under two scenarios, considering i) modelled atmospheric data and ii) atmospheric measurements. The error estimation was performed first by means of the RMSE and MBE estimators calculated for two-months of temperature and AH modelled HYSPLIT radiosondes and MWR. The aim was to quantify the error propagation on the two atmospheric variables involved in the RCS w

(

  DCN). In this work, DCN measurements were carried out under day-time and night-time conditions using the termination hood accessory delivered with the ceilometer shipping for covering the instrument to avoid external light contamination, however some reflections in the near field are still remained (below 50m). The DCN was regularly measured twice per week during two weeks, taking samples of 30 minutes at day and night time.During the analyzed period (October-November 2019) DCN does not show a high variability between days, therefore we focused our attention on the day-time and night-time analysis. As example, Fig.2apresents the day-time(gray line) and night-time (black line) mean DCN profiles measured from 0.010 km to 15.4 km agl for the 11 th October 2019. The bias between the profiles (not shown) in the first kilometers after full overlap (from 0.80m to 1km agl) is close to zero, but it increased rapidly with height: bias from 1 to 2 km agl increased from 2.5 to 5.0 a.u, from 2 to 6 km agl increased from 5 to 25 a.u and above this height the bias reached up values above 300 a.u. The shape

Figure 1 :

 1 Figure 1: Block diagram of the methodology

Figure 2 :

 2 Figure 2: Example measurements corresponding to 11 th October 2019. Here, the ceilometer DCN correction is shown. a) The DCN measurements are presented in black: night-time and gray: day-time; b) the smoothed DCN smoothed profiles (for the sake of clarity, but not considered in the calculations); c) the RCS raw signal 1h-averaged, and d) the RCS raw DCN corrected (smoothed)

Figure 3 :

 3 Figure 3: The figure presents the noise pre-processing applied to obtain the RCS * =RSC raw -(DC+BG). The measurements were took on 11 th October 2019. In the upper part are presented the results at 7 UTC and the bottom at 10UTC. The four BG ranges are shown. In green line is related to the correction for BG1, orange to BG2, blue to BG3 and cyan to BG4. The standard deviation calculated corresponds to 500m-average in order to gain visibility at each step of the signal pre-processing.

Figure 4 :

 4 Figure 4: Water vapor correction corresponding to 11 th October 2019 at 06:00 UTC. a) Blue line represents the AH profile, b) RCS w with/without water vapor correction, red and gray line respectively are shown , and c) is the bias between RCS * and RCS w

(

  LR m ) from theory. The other terms like particle lidar ratio (LR p ), and the reference height (r ref ) are calculated during the different steps involved in the modified[START_REF] Cazorla | Near-real-time processing of a ceilometer network assisted with sunphotometer data: monitoring a dust outbreak over the iberian peninsula[END_REF] algorithm to Vaisala ceilometer.An iterative Klett inversion process is performed for 240 values of LR ranging from 35 to 150 Sr each 0.5 sr, and for all possible height reference points inside Rayleigh zone found above, in order to retrieve β p and α p . The goal with this iterative process is to run Klett until find a combination of LR and reference Rayleigh height that makes comparable the integrated α p profile with the interpolated AOD 910 from sun photometer at the closest UTC time. The retrievals are considered as successful when the difference between integrated α p profile and AOD 910 is lower than 0.001.

Figure 5 :

 5 Figure 5: Example case of Klett inversion performed on 11 th October 2019 for the time-frames that sun photometer has LR data availability. The AOD 910 was interpolated from AERONET data at the same UTC time of the 1h-averaged RCS w . The black dotted line in a) refers to the Rayleigh fit profile and in colors are presented in a) the RCS w profiles used as input for the retrieval, while in b) is represented the α p resulting from the inversion

  Figure 6a,b shows in red line the temperature profiles of RMSE and MBE calculations, pointing out that below 2 km agl the RMSE reach the lowest values (lower than 4 • C),

Figure 6 :

 6 Figure 6: It is presented the error calculation in terms of RMSE and MBE, panel a) and b) respectively. For both panels, red line refers to temperature and blue line to AH.

Figure 7 :

 7 Figure 7: The panel presents the different steps in the example case on 11 th October 2019 at 8 UTC for evaluation the relative error (RE) in the inversion process. a) In red lines is represented the temperature and blue line refers to AH, b) is the error committed in the RCS w signal, and c) RE of the β p (red line), and α p (black line)

Table 1 :

 1 Klett inversions on 11 th October 2019.

	UTC Klett AOD Klett LR (sr)
	8	0.039	72.2
	9	0.045	58.1
	10	0.052	68.7
	12	0.059	83.6
	13	0.060	99.0
	15	0.061	87.9

2019 at 06:00 UTC and cases 5 and 6 to 21 st November 2019 at 09:00 and 12:00 UTC. The averages were calculated for 3 different layers for center the attention in the aerosol products: from 0 to 0.5 km agl, from 0.5 to 1 km agl, from 1 to 2 km agl. The inversion products (β p and α p ) were analyzed until 2 km agl.

Table2 shows that mean temperature RE were always below 4 % with lower standard deviation (SD). However, in the first 0.5km agl case 2 and 3 presented mean values up to 4 %. On the contrary, AH presented higher mean error values mainly above 3 km agl as it was expected from the analysis performed in the previous section with mean errors higher than 40 % (not shown here). Centering the attention in the first two kilometres of the atmosphere, the errors in AH were below 30 %, except for case 5,6 with mean errors up to 55 %. As it was seen on Sec.3, the AH is deeply linked with the transmittance term at correcting the RCS * signal, therefore the relatively large differences in AH profiles are contributing to RCS w profiles. According to our results, mean RE up to 30 % in AH might cause errors in RCS w that reach 10 % in the first 2 km agl, and like the amount of water vapor is highly variable in the first two kilometers of the atmosphere, nearly-real measurements of AH will improve the ceilometer signal significantly, otherwise calculations without using real ambient measurements might lead to considerable error increase.

The retrieved products β p and α p have a RE that is influenced firstly by the input profiles which are Rayleigh and RCS w (including the noises remained), and secondly the accurate atmospheric reference to start the backward inversion, which is linked with the sunphotometric data for setting the correct LR, among other errors associated to the algorithm calculation itself. From Table2, it is possible to see that in general, the mean RE is lower than 3 % below 2 km agl for β p , meanwhile the errors in α p can reach 25 % within this height. Some mean errors for α p in cases 4,5,6 are quite larger than other cases, and this might be linked with the relative elevated AH errors (see on Table2). Considering this error propagation in α p , we found that without using co-located atmospheric measurements, this inversion product can be also estimated but considering an addition up to 25 % of RE for Vaisala CL51 ceilometers, avoiding those cases where AH RE was extremely high. The slight increase in the standard deviation errors reported for the RCS w in the last air volume evaluated from 1 to 2 km agl (Table2) are associated to the correction. The errors were calculated for 488 samples covering different atmospheric scenarios, showing that for temperatures below 2km agl the RMSE is lower than 4 • C, and above is up to 6 • C, whereas AH presented an error up to 2 gm -3 , which in terms of water vapor correction is significant. The MBE shown that temperature modelled profiles were always underestimating MWR ones (up to -1 • C), and AH MBE was high variable, underestimating the MWR measures from 1 km to 2.3 km agl and from 5.8 km agl to the end of the profile (reaching up -0.5 gm -3 ) and the overestimation reached up 1.5 gm -3 from ground to 0.2 km agl and from 3 km to 5.8 km agl.

In summary, we found that the use of modelled data instead atmospheric measurements is primarily influencing the RCS w , and then the inversion products obtained from the ceilometer. As a result of that, the percentual relative error estimation within the critical zone for aerosol inversion (first 2km agl), the temperature presented errors below 4%, while in the first 2 km agl AH were below 25 %. These errors affected the RCS w signal driving to an error up to 9 % in the first 2km agl. The errors propagated on β p during Klett calculation were lower than 2.2 %, leading to an error in α p up to 25 %. Therefore, one can conclude that the use of atmospheric modelled data instead of measurements for water vapor correction on vaisala CL51 ceilometers will lead large errors on inversion products i.e. α p .

Finally, the Klett algorithm could be improved by determining the calibration constant of the instrument tackling one of the E-PROFILE objectives. This work can be developed in further studies by having a larger ceilometer database for making a seasonal analysis of the calibration constant as function of the internal temperature in a semi-automatic way at ONERA site. In addition, the knowledge of the full overlap height of the system could also help us to improve the inversion products, therefore a further measurement campaign with ground based lidar operating co-located to the ceilometer might be an ideal solution to have the ceilometer fully characterized.


Xuan Son Nguyen

Luc Brun

Olivier Lézoray

Sébastien Bougleux

Learning Recurrent High-order Statistics for Skeleton-based Hand Gesture Recognition

High-order statistics have been proven useful in the framework of Convolutional Neural Networks (CNN) for a variety of computer vision tasks. In this paper, we propose to exploit high-order statistics in the framework of Recurrent Neural Networks (RNN) for skeleton-based hand gesture recognition. Our method is based on the Statistical Recurrent Units (SRU), an un-gated architecture that has been introduced as an alternative model for Long-Short Term Memory (LSTM) and Gate Recurrent Unit (GRU). The SRU captures sequential information by generating recurrent statistics that depend on a context of previously seen data and by computing moving averages at different scales. The integration of high-order statistics in the SRU significantly improves the performance of the original one, resulting in a model that is competitive to state-of-the-art methods on the Dynamic Hand Gesture (DHG) dataset, and outperforms them on the First-Person Hand Action (FPHA) dataset.

I. INTRODUCTION

Skeleton-based hand gesture recognition has been an active research topic in recent years with many potential applications in assisted living, human-robot interaction, sign language interpretation, smart home control interface. In this work, we focus on deep learning techniques [START_REF] Liu | Recognizing Human Actions as The Evolution of Pose Estimation Maps[END_REF], [START_REF] Wang | Action Recognition Based on Joint Trajectory Maps Using Convolutional Neural Networks[END_REF], [START_REF] Du | Hierarchical Recurrent Neural Network for Skeleton Based Action Recognition[END_REF], [START_REF] Wang | Modeling Temporal Dynamics and Spatial Configurations of Actions Using Two-Stream Recurrent Neural Networks[END_REF], [START_REF] Liu | Global Context-Aware Attention LSTM Networks for 3D Action Recognition[END_REF], [START_REF] Shahroudy | NTU RGB+D: A Large Scale Dataset for 3D Human Activity Analysis[END_REF] for solving the task as they have shown superior performance over their counterparts based on hand-crafted features [START_REF] Luo | Group Sparsity and Geometry Constrained Dictionary Learning for Action Recognition from Depth Maps[END_REF], [START_REF] Yang | EigenJoints-based Action Recognition Using Naive-Bayes-Nearest-Neighbor[END_REF], [START_REF] Vemulapalli | Human Action Recognition by Representing 3D Skeletons as Points in a Lie Group[END_REF]. Recently, modeling feature distribution with highorder statistics in deep neural networks has been proven effective. While this idea has been proposed in a number of CNN-based approaches [START_REF] Lin | Bilinear CNN Models for Fine-Grained Visual Recognition[END_REF], [START_REF] Li | Is Second-order Information Helpful for Large-scale Visual Recognition?[END_REF], [START_REF] Wang | G2DeNet: Global Gaussian Distribution Embedding Network and Its Application to Visual Recognition[END_REF], very few RNN-based approaches [START_REF] Chakraborty | A Statistical Recurrent Model on the Manifold of Symmetric Positive Definite Matrices[END_REF] exploit statistical information higher than first-order. One of the main challenges with such approaches is that the use of high-order statistics results in data lying on a manifold. This requires carefully-designed operations taking into account the geometric structure of data processed by the network [START_REF] Huang | Building Deep Networks on Grassmann Manifolds[END_REF].

In this work, we propose an approach that combines the ideas of the SRU [START_REF] Oliva | The Statistical Recurrent Unit[END_REF] for learning long term information in sequences and ST-TS-HGR-NET [START_REF] Nguyen | A Neural Network Based on SPD Manifold Learning for Skeleton-based Hand Gesture Recognition[END_REF], recently introduced for hand gesture recognition. In particular, we introduce highorder statistics into the SRU for improving hand gesture recognition. The motivation to use the SRU is that it has a simple un-gated architecture while being competitive to more sophisticated LSTM and GRU alternatives [START_REF] Oliva | The Statistical Recurrent Unit[END_REF]. This is particularly useful when we want to introduce high-order statistics into the model, since it will simplify the design of operations taking into account the geometric structure of the resulting data. Note that SRU has not been proposed for hand gesture recognition in previous works. Our model is designed to capture the temporal dependency of matrices encoding firstorder and second-order statistics of hand joints' coordinates. Moreover, we consider learning statistics resulting from the "covariance" of those matrices that significantly improves the performance. As a result, our method is superior to SRU on the DHG and FPHA datasets. Compared to ST-TS-HGR-NET, our method gives competitive results on the DHG dataset and outperform it on the FPHA dataset while using far less parameters.

II. RELATED WORKS

A. Skeleton-Based Gesture Recognition

Existing approaches can be roughly classified into handcrafted feature based approaches or deep learning approaches. Hand-crafted feature based approaches compute joint features from the relative position between pairs of joints [START_REF] Luo | Group Sparsity and Geometry Constrained Dictionary Learning for Action Recognition from Depth Maps[END_REF], [START_REF] Smedt | Skeleton-Based Dynamic Hand Gesture Recognition[END_REF], [START_REF] Wang | Mining Actionlet Ensemble for Action Recognition with Depth Cameras[END_REF], [START_REF] Yang | EigenJoints-based Action Recognition Using Naive-Bayes-Nearest-Neighbor[END_REF], the relative position of 4-tuples of joints [START_REF] Evangelidis | Skeletal Quads: Human Action Recognition Using Joint Quadruples[END_REF], Gram matrices [START_REF] Zhang | Efficient Temporal Sequence Comparison and Classification Using Gram Matrix Embeddings on a Riemannian Manifold[END_REF], or the rotations and translations describing the 3D geometric relationships between body parts [START_REF] Vemulapalli | Human Action Recognition by Representing 3D Skeletons as Points in a Lie Group[END_REF]. Deep learning approaches based on CNN [START_REF] Devineau | Deep Learning for Hand Gesture Recognition on Skeletal Data[END_REF], [START_REF] Ke | A New Representation of Skeleton Sequences for 3D Action Recognition[END_REF], [START_REF] Liu | Enhanced Skeleton Visualization for View Invariant Human Action Recognition[END_REF], [START_REF] Liu | Recognizing Human Actions as The Evolution of Pose Estimation Maps[END_REF], [START_REF] Nez | Convolutional Neural Networks and Long Short-Term Memory for Skeleton-based Human Activity and Hand Gesture Recognition[END_REF], [START_REF] Nguyen | A Neural Network Based on SPD Manifold Learning for Skeleton-based Hand Gesture Recognition[END_REF], [START_REF] Wang | Action Recognition Based on Joint Trajectory Maps Using Convolutional Neural Networks[END_REF], RNN [START_REF] Du | Hierarchical Recurrent Neural Network for Skeleton Based Action Recognition[END_REF], [START_REF] Wang | Modeling Temporal Dynamics and Spatial Configurations of Actions Using Two-Stream Recurrent Neural Networks[END_REF] and LSTM [START_REF] Liu | Spatio-Temporal LSTM with Trust Gates for 3D Human Action Recognition[END_REF], [START_REF] Liu | Global Context-Aware Attention LSTM Networks for 3D Action Recognition[END_REF], [START_REF] Nez | Convolutional Neural Networks and Long Short-Term Memory for Skeleton-based Human Activity and Hand Gesture Recognition[END_REF], [START_REF] Shahroudy | NTU RGB+D: A Large Scale Dataset for 3D Human Activity Analysis[END_REF], [START_REF] Weng | Deformable Pose Traversal Convolution for 3D Action and Gesture Recognition[END_REF] have demonstrated impressive results. Some works have shown the interest of deep learning on manifolds for action and hand gesture recognition, e.g., SPD manifolds [START_REF] Huang | A Riemannian Network for SPD Matrix Learning[END_REF], [START_REF] Nguyen | A Neural Network Based on SPD Manifold Learning for Skeleton-based Hand Gesture Recognition[END_REF], Grassmann manifolds [START_REF] Huang | Building Deep Networks on Grassmann Manifolds[END_REF], and Lie groups [START_REF] Huang | Deep Learning on Lie Groups for Skeleton-Based Action Recognition[END_REF]. Since a hand or a body skeleton can be naturally represented as a graph (see e.g., Fig. 1), deep learning on graphs has also been applied to skeleton-based action recognition [START_REF] Li | Spatio-Temporal Graph Convolution for Skeleton Based Action Recognition[END_REF], [START_REF] Yan | Spatial Temporal Graph Convolutional Networks for Skeleton-Based Action Recognition[END_REF].

B. Deep Learning with Second-order Statistics

Modeling feature distribution with high-order statistics has been proven effective in hand-crafted feature based approaches for various vision tasks. Inspired by this idea, some recent CNNs [START_REF] Ionescu | Matrix Backpropagation for Deep Networks with Structured Layers[END_REF], [START_REF] Li | Is Second-order Information Helpful for Large-scale Visual Recognition?[END_REF], [START_REF] Lin | Bilinear CNN Models for Fine-Grained Visual Recognition[END_REF], [START_REF] Yu | Statistically-motivated Second-order Pooling[END_REF] generate image representations by performing second-order pooling (as a covariance matrix) after the last convolutional layers of the networks. They have shown to significantly outperform classical networks [START_REF] He | Deep Residual Learning for Image Recognition[END_REF], [START_REF] Huang | Densely Connected Convolutional Networks[END_REF], [START_REF] Szegedy | Going Deeper with Convolutions[END_REF] based on first-order statistics (mean vector) for producing image representations. Pooling strategies combining first-order and second-order statistics [START_REF] Dai | FASON: First and Second Order Information Fusion Network for Texture Recognition[END_REF], [START_REF] Nguyen | A Neural Network Based on SPD Manifold Learning for Skeleton-based Hand Gesture Recognition[END_REF], [START_REF] Wang | G2DeNet: Global Gaussian Distribution Embedding Network and Its Application to Visual Recognition[END_REF] have also demonstrated better performance than those based only on first-order or second-order statistics. Wang et al. [START_REF] Wang | Global Gated Mixture of Second-order Pooling for Improving Deep Convolutional Neural Networks[END_REF] proposed a second-order pooling method capable of capturing complex and non-unimodal feature distributions. This relaxes the unimodal distribution assumption required by the above works and thus helps increasing the representation ability of CNNs. Kernel-based methods [START_REF] Cai | Higher-Order Integration of Hierarchical Convolutional Activations for Fine-Grained Visual Categorization[END_REF], [START_REF] Cui | Kernel Pooling for Convolutional Neural Networks[END_REF], [START_REF] Gao | Compact Bilinear Pooling[END_REF] reveal the link between high-order pooling and kernel machines, achieving compact covariance matrix-based representations. While most approaches introduced second-order pooling layers operating at the end of CNNs, Gao et al. [START_REF] Gao | Global Second-order Pooling Convolutional Networks[END_REF] proposed a second-order pooling block that can be inserted after any convolutional layer. This allows CNNs to capture second-order statistics in intermediate layers. Similarly, the works of [START_REF] Li | Factorized Bilinear Models for Image Recognition[END_REF], [START_REF] Wang | SORT: Second-Order Response Transform for Visual Recognition[END_REF], [START_REF] Zoumpourlis | Non-linear Convolution Filters for CNN-Based Learning[END_REF] focus on modeling feature interaction in intermediate layers of CNNs. Some works [START_REF] Chakraborty | A Statistical Recurrent Model on the Manifold of Symmetric Positive Definite Matrices[END_REF], [START_REF] Dong | Deep Manifold Learning of Symmetric Positive Definite Matrices with Application to Face Recognition[END_REF], [START_REF] Engin | DeepKSPD: Learning kernelmatrix-based spd representation for fine-grained image recognition[END_REF], [START_REF] Huang | A Riemannian Network for SPD Matrix Learning[END_REF] aim at learning covariance matrices that preserve the same geometric structure of the input data.

III. STATISTICAL RECURRENT UNIT (SRU)

Oliva et al. [START_REF] Oliva | The Statistical Recurrent Unit[END_REF] proposed an un-gated unit referred to as the statistical recurrent unit (SRU) for learning long term information in sequences. The SRU captures sequential information by generating recurrent statistics that depend on a context of previously seen data and by computing moving averages at different scales. Let x 1 , x 2 , . . . , x T be an input sequence on R n . The update rules for the SRU are as follows:

R t = ReLU(W r χ χ χ t-1 + b r ) (1)
P t = ReLU(W p R t + W x x t + b p ) (2)
∀α ∈ A, χ χ χ

(α) t = αχ χ χ (α) t-1 + (1 -α)P t (3)
O t = ReLU(W o χ χ χ t + b o ) (4)
where } α∈A and then used to create an output of the network at frame t.

Note that the moving averages are dependent not only on the current input but also on the statistics computed from values of previous points. This approach is different from treating the sequence as a set of i.i.d. points drawn from some distribution and marginalizing out time, which loses temporal information. While the SRU is based on an un-gated architecture, it has been shown to be competitive with more sophisticated LSTM and GRU alternatives [START_REF] Oliva | The Statistical Recurrent Unit[END_REF].

IV. PROPOSED APPROACH

We propose in Section IV-A a method for constructing our network input which takes into account the correlations of neighboring hand joints. In Section IV-B, a sequence model based on SRU and high-order statistics is then introduced. Based on the proposed model, we show how to build a network architecture for hand gesture recognition in Section IV-C. Finally, an interpretation of the model from a statistical point of view is provided in Section IV-D.

A. Network Input

In this work, we use skeletal data obtained by an Intel RealSense camera for hand gesture recognition. The hand joints' positions at each frame are illustrated in Fig. 1. Given a sequence of hand joints' 3D coordinates, we first remove the two joints at the palm and wrist (if any) and then renumber the joint indices from 1 to N (N = 20). The construction of our network input at a given frame t is illustrated in Fig. 2. Fig. 2(a) represents the 3-dimensional array containing the hand joints' 3D coordinates at frame t where each number indicates a joint index. Each joint i, i = 1, . . . , N belonging to a finger has at most two neighboring joints belonging to the same finger, referred to as the lower and upper neighbors of i. For a hand joint with no lower neighbor [START_REF] Liu | Recognizing Human Actions as The Evolution of Pose Estimation Maps[END_REF][START_REF] Liu | Global Context-Aware Attention LSTM Networks for 3D Action Recognition[END_REF][START_REF] Vemulapalli | Human Action Recognition by Representing 3D Skeletons as Points in a Lie Group[END_REF][START_REF] Chakraborty | A Statistical Recurrent Model on the Manifold of Symmetric Positive Definite Matrices[END_REF][START_REF] Smedt | Skeleton-Based Dynamic Hand Gesture Recognition[END_REF] or no upper neighbor [START_REF] Wang | Modeling Temporal Dynamics and Spatial Configurations of Actions Using Two-Stream Recurrent Neural Networks[END_REF][START_REF] Yang | EigenJoints-based Action Recognition Using Naive-Bayes-Nearest-Neighbor[END_REF][START_REF] Wang | G2DeNet: Global Gaussian Distribution Embedding Network and Its Application to Visual Recognition[END_REF][START_REF] Nguyen | A Neural Network Based on SPD Manifold Learning for Skeleton-based Hand Gesture Recognition[END_REF][START_REF] Zhang | Efficient Temporal Sequence Comparison and Classification Using Gram Matrix Embeddings on a Riemannian Manifold[END_REF], we duplicate its unique neighbor so that it has two neighbors. The 3-dimensional array in Fig. 2(b) is constructed by replacing the coordinates of each joint in Fig. 2(a) by those of its lower neighbor. Similarly, the 3-dimensional array in Fig. 2(c) is constructed by replacing the coordinates of each joint in Fig. 2(a) by those of its upper neighbor. Our network input at frame t is then formed by concatenating the three 3-dimensional arrays along the depth dimension (Fig. 2(d)). Thus, the feature vector of a hand joint has 9 dimensions. This method for constructing the network input allows to take into account the correlations of neighboring hand joints and has been shown to greatly improve recognition accuracy (see Section V).

B. SRU Based on High-order Statistics (SRU-HOS)

Differently from previous approaches that are based only on first-order and second-order statistics, we propose to combine those statistics with higher-order statistics computed from

R f,k t P f,k t X f,k t X f,k t T k +1 χ χ χ f,k t 1 O f,k t R t P t χ χ χ t 1 O t X t T k +1:t
T 2 = t 1 + t 2 -1.
them. For each finger f = 1, . . . , 5, we consider two statistics (χ χ χ f,k t ) k=1,2 at frame t that are updated using the following equations (see Fig. 3):

R f,k t = ReEig α∈A (w k,α r ) 2 α∈A (w k,α r ) 2 χ χ χ f,k,α t-1 (5)
P f,k t = ReEig (w k p ) 2 (w k p ) 2 + (w k x ) 2 R f,k t + (w k x ) 2 (w k p ) 2 + (w k x ) 2 h k (X f t ) (6)
∀α ∈ A, χ χ χ f,k,α t = αχ χ χ f,k,α t-1 + (1 -α)P f,k t (7) O f,k t = ReEig α∈A (w k,α o ) 2 α∈A (w k,α o ) 2 χ χ χ f,k,α t , ( 8
)
where

A is a set of different scales, w k,α r , w k p , w k x , w k,α o ∈ R are the parameters of the model, X f t ∈ R 9×|J f |
is the matrix at frame t representing the 3D coordinates of hand joints belonging to finger f (Section IV-A), J f is the set of hand joints belonging to finger f , χ χ χ f,k,α 0 are the initial states, h 1 (X f t ) and h 2 (X f t ) compute statistics from frames tt 1 + 1 to t and frames tt 1t 2 + 2 to t, t 1 , t 2 > 0 are two constants, and ReEig(Y) = U max( I, V)U T , where Y is a symmetric positive definite (SPD) matrix and Y = UVU T its eigen-decomposition, is a rectification threshold, I is the identity matrix, max( I, V) is a diagonal matrix whose diagonal elements are defined as:

(max( I, V))(i, i) = V(i, i) if V(i, i) > if V(i, i) ≤ . (9)
The function ReEig was proposed in [START_REF] Huang | A Riemannian Network for SPD Matrix Learning[END_REF] for introducing non-linear transformation of SPD matrices. Let X t = [x 1 t , . . . , x N t ], where x i t ∈ R 9 represents the feature vector of joint i at frame t, i = 1, . . . , N . Assuming that x i j , i ∈ J f , j = tt 1 + 1, . . . , t are independent and identically distributed samples from the following Gaussian distribution:

N (x; µ µ µ f,1 t , Σ Σ Σ f,1 t ) = 1 |2πΣ Σ Σ f,1 t | 1 2 exp - 1 2 (x -µ µ µ f,1 t ) T (Σ Σ Σ f,1 t ) -1 (x -µ µ µ f,1 t ) , (10)
where |.| is the determinant, µ µ µ f,1 t is the mean vector and Σ Σ Σ f,1 t is the covariance matrix:

µ µ µ f,1 t = 1 |J f |t 1 t j=t-t1+1 i∈J f x i j , Σ Σ Σ f,1 t = 1 |J f |t 1 t j=t-t1+1 i∈J f (x i j -µ µ µ f,1 t )(x i j -µ µ µ f,1 t ) T . Then h 1 (X f t ) is designed to capture the Gaussian distribution N (x; µ µ µ f,1 t , Σ Σ Σ f,1 t ) [47]: h 1 (X f t ) = Σ Σ Σ f,1 t + µ µ µ f,1 t (µ µ µ f,1 t ) T µ µ µ f,1 t (µ µ µ f,1 t ) T 1 . (11)
By introducing h 1 (X f t ) in Eq. ( 6), the proposed model learns the first-order (mean) and second-order (covariance) statistics of the hand joints' feature vector. To strengthen the learning capabilities of the model, it also learns higherorder statistic h 2 (X f t ) defined by the "covariance" of h 1 (X f t ). Since h 1 (X f t ) lies on the manifold of SPD matrices, it is projected onto its tangent space before the covariance can be measured. The projection is performed by the matrix function log p (Y) = U log(V)U T [START_REF] Arsigny | Geometric Means in a Novel Vector Space Structure on Symmetric Positive Definite Matrices[END_REF], where Y is a SPD matrix and Y = UVU T is its eigen-decomposition. Then the statistic h 2 (X f t ) is defined as:

µ µ µ f,2 t = 1 t2 t i=t-t2+1 vl(h 1 (X f i )), (12)
h 2 (X f t ) = 1 t2 t i=t-t2+1 (vl(h 1 (X f i )) -µ µ µ f,2 t ) • (vl(h 1 (X f i )) -µ µ µ f,2 t ) T , (13)
where vl(•) = vec(log(•)), vec : R d×d → R d(d+1)/2 the special vectorization operator [START_REF] Tuzel | Pedestrian Detection via Classification on Riemannian Manifolds[END_REF] that transforms any symmetric matrix

Y ∈ R d×d into a vector vec(Y) = [Y(u, u), ( √ 2Y(u, v)) v=u+1,...,d ] T u=1,...,d .
Let M be the index of the last frame of the sequence, X f be the matrix representing the 3D coordinates of hand joints belonging to finger f , O f,k M , k = 1, 2 be the final outputs of our model. We can then write the final outputs of our model as a function of the input and the model parameters as follows:

{O f,k M } k=1,2 =SRU-HOS(X f , {χ χ χ f,k,α 0 } k=1,2,α∈A , {w k p } k=1,2 , {w k x } k=1,2 , {w k,α r } k=1,2,α∈A , {w k,α o } k=1,2,α∈A ) (14)
Õs,f,2

Ms

Õs,f,1

Ms O s,f,2 t-T2+1 O s,f,1 t-T2+1 O s,f,1 t-1 O s,f,2 t-1 O s,f,2 Ms O s,f,1 Ms O s,f,2 t O s,f,1 t X t-T2+1 X t-1 X Ms X t vl(•)

C. The Proposed Network

Our proposed network is illustrated in Fig. 4. To better capture temporal dependencies of a skeleton sequence, we create six sub-sequences by dividing the original sequence into one sequence, then two sequences of equal length, and three sequences of equal length. For each sub-sequence s, s = 1, . . . , 6, and a finger f, f = 1, . . . , 5, we create a model for learning statistics from s associated with f . Denote by X s,f the matrix representing the 3D coordinates of hand joints belonging to finger within sub-sequence s, χ χ χ s,f,k,α 0 are the initial states, w s,k,α r , w s,k p , w s,k x , w s,k,α o ∈ R are the parameters of the model. Then the final outputs of the model can be given by:

{O s,f,k Ms } k=1,2 =SRU-HOS(X s,f , {χ χ χ s,f,k,α 0 } k=1,2,α∈A , {w s,k p } k=1,2 , {w s,k x } k=1,2 , {w s,k,α r } k=1,2,α∈A , {w s,k,α o } k=1,2,α∈A ) (15)
For a more accurate classification, the outputs

O s,f,1 Ms ∈ R 10×10 and O s,f,2 Ms ∈ R 55×55 (Eq. (8)) are tuned w.r.t. parameters W s,f,1 ∈ R 10×10 and W s,f,2 ∈ R 55×55 : Õs,f,1 Ms = W s,f,1 O s,f,1 Ms (W s,f,1 ) T , (16)
Õs,f,2

Ms = W s,f,2 O s,f,2 Ms (W s,f,2 ) T . ( 17
)
Both matrices W s,f,k , k = 1, 2 have full row-rank so that Õs,f,k

Ms

is SPD as O s,f,k Ms . These weights can be optimized using a variant of stochastic gradient descent on Stiefel manifolds [START_REF] Huang | A Riemannian Network for SPD Matrix Learning[END_REF]. Since O s,f,2 Ms encodes the "covariance" of h 1 (X s,f t ), these outputs are grouped into the matrix Y s,f to obtain the final representation of sub-sequence s associated with finger f as follows:

Y s,f = Õs,f,2 Ms + vl( Õs,f,1 Ms )vl( Õs,f,1 Ms ) T vl( Õs,f,1 Ms ) vl( Õs,f,1 Ms ) T 1 , (18)
We transform Y s,f into an Euclidean space and then vectorize the resulting matrix using the vl(.) function to create a representation of sub-sequence s associated with finger f . A global representation of the original sequence is obtained by concatenating all resulting vectors for s = 1, . . . , 6 and f = 1, . . . , 5, i.e. [vl(Y 1,1 ) T , . . . , vl(Y 6,5 ) T ] T . Finally, the network is trained by passing this vector to a fully connected layer and by minimizing cross-entropy (see Fig. 4).

It is worth pointing out the differences between our network and ST-TS-HGR-NET. First, both networks capture weak orders of frames by computing statistics from six temporal subsequences and then combining them to represent the original sequence. However, each sub-sequence in ST-TS-HGR-NET is represented by the mean and covariance computed from statistics of all frames of the sub-sequence without taking into account the temporal order of frames. In contrast, this order is used in Eq. ( 6) of our network where P s,f,k t is computed from the current input X s,f t and the statistics of the previous data R s,f,k t . Thus, temporal information could be better preserved by our network than ST-TS-HGR-NET. Second, in our network, different sets of weights are learned for different sub-sequences but these weights are shared for different fingers. In ST-TS-HGR-NET, a set of weights is learned for each sub-sequence associated with each finger. This results in a large set of parameters to be trained compared to our network (see Section V-C).

D. Interpretation of the Model

From a statistical point of view, our model can be interpreted as an extension of the SRU model. Let (y t ) t=1,...,M be an input sequence of real-valued points and let (γ t ) t=1,...,M be the statistics on these points. In a SRU, the statistic γ t on the t-th point y t is defined as a function of y t and the statistic γ t-1 on the previous point: γ t = γ(y t , γ t-1 ) with γ 0 an initial constant vector. In our model, two statistics γ 1 t and γ 2 t are computed for t = t 1 + t 2 -1, . . . , M . Each statistic is defined as a function of the statistic on the previous point, and the statistic generated by a function h or g from y t and previous points:

γ 1 t = γ(h(y t-t1+1 , . . . , y t ), γ 1 t-1 ), (19)
γ 2 t = γ(g(h(y t-t1-t2+2 , . . . , y t-t2+1 ), . . . , h(y t-t1+1 , . . . , y t )), γ 2 t-1 ).

(

) 20
As with an SRU, the statistics produced by Eq. ( 19) can perfectly encode a sequence if our model is given enough dimensions. For example, consider a sequence y 1 , . . . , y M ∈ R + and statistics γ t = (0, . . . , 0) for t = 1, . . . , t 1 -1 and γ t = (0, . . . , M t1 y t-t1+1 , . . . , M t1 y t , 0, . . .) for t = t 1 , . . . , M . The values at indices t-t 1 +1, . . . , t are M t1 y t-t1+1 , . . . , M t1 y t , respectively. The average of these statistics is given by:

The average can be used to recover the original sequence. Since the statistics generated by h depend on the current and previous points, the number of samples for the estimate of Σ Σ Σ s,f,1 t is increased. It is particularly useful in our method to avoid ill-conditioned matrices (Eq. ( 11)).

V. EXPERIMENTS

Our network, referred to as SRU-HOS-NET, is validated on the Dynamic Hand Gesture (DHG) dataset [START_REF] Smedt | Skeleton-Based Dynamic Hand Gesture Recognition[END_REF], [START_REF] Smedt | 3D Hand Gesture Recognition Using a Depth and Skeletal Dataset[END_REF] and the First-Person Hand Action (FPHA) dataset [START_REF] Garcia-Hernando | First-Person Hand Action Benchmark with RGB-D Videos and 3D Hand Pose Annotations[END_REF]. The batch size and the learning rate were set to 40 and 0.9, respectively. The set of values for α was set similarly to [START_REF] Oliva | The Statistical Recurrent Unit[END_REF], i.e., A = {0.01, 0.25, 0.5, 0.9, 0.99}. The rectification threshold for the ReEig(•) function was set to 10 -4 [START_REF] Huang | A Riemannian Network for SPD Matrix Learning[END_REF]. For sequences having more or less than 150 frames, we used interpolation in order to fix the number of frames of all sequences to 150. Gesture recognition was performed by training a SVM classifier using the output of SRU-HOS-NET (the vector obtained before the fully connected layer) at epoch 50 (see Fig. 4). In our experiments, we found that this method for classifying gestures gave competitive results compared to those obtained by training the network with a softmax layer. Indeed, if one uses a classical FC/SoftMax Layer for the classification, the networks usually converges in more than 200 epochs. However similar results can be obtained with a SVM taking the output of the network trained during four times less epochs. We used the LIBLINEAR library [START_REF] Fan | LIBLINEAR: A Library for Large Linear Classification[END_REF] to train our SVM classifier with L2-regularized L2-loss (dual) and default parameter settings.

A. Datasets and Experimental Protocols

DHG dataset. It contains 14 gestures performed in two ways: using one finger and the whole hand. Each gesture is executed several times by different actors. The gestures are subdivided into categories of fine and coarse : grab (fine), tap (coarse), expand (fine), pinch (fine), rotation clockwise (rotcw, fine), rotation counterclockwise (rot-ccw, fine), swipe right (swipe-r, coarse), swipe left (swipe-l, coarse), swipe up (swipeu, coarse), swipe down (swipe-d, coarse), swipe x (swipe-x, coarse), swipe v (swipe-v, coarse), swipe + (swipe-+, coarse), shake (coarse). The dataset provides the 3D coordinates of 22 hand joints as illustrated in Fig. 1. It has been split into 1960 train sequences (70% of the dataset) and 840 test sequences (30% of the dataset) [START_REF] Smedt | 3D Hand Gesture Recognition Using a Depth and Skeletal Dataset[END_REF]. pour juice bottle, pour liquid soap, pour milk, pour wine, prick, put salt, put sugar, put tea bag, read letter, receive coin, scoop spoon, scratch sponge, sprinkle, squeeze paper, squeeze sponge, stir, take letter from envelope, tear paper, toast wine, unfold glasses, use calculator, use flash, wash sponge, write. Action sequences present high inter-subject and intra-subject variability of style, speed, scale, and viewpoint. The dataset provides the 3D coordinates of 21 hand joints as DHG dataset except for the palm joint. We used the 1:1 setting proposed in [START_REF] Garcia-Hernando | First-Person Hand Action Benchmark with RGB-D Videos and 3D Hand Pose Annotations[END_REF] with 600 action sequences for training and 575 for testing.

B. Ablation Study

In this section, we investigate the impact of different components on our network's performance. We set the default values of t 2 to 15. Time windows t 1 . Tab. I (column t 1 ) shows the impact of the time window t 1 on SRU-HOS-NET, where t 1 is set to t 1 = 3, 10, 20. As can be observed, SRU-HOS-NET achieves the best results with t 1 = 20 and t 1 = 3 on DHG and FPHA datasets, respectively. Results show that combining the current point with a sufficient number of previous points for estimating h1 (•) can lead to non-trivial performance gain. In the following, we set t 1 to 20 and 3 for experiments on DHG and FPHA datasets, respectively.

Statistics

Time windows t 2 . The impact of the time window t 2 on SRU-HOS-NET is shown in Tab. I (column t 2 ). Here we test with three different settings of t 2 : t 2 = 10, 15, 20. For both the datasets, SRU-HOS-NET gives the best results with t 2 = 15.

Network input. The impact of the feature concatenation method (section IV-A) can be seen in Tab. I (column Feature concatenation), where the performance of SRU-HOS-NET drops significantly when feature concatenation is not used, i.e., the hand joints' coordinates are fed directly as input of SRU-HOS-NET. Results indicate that feature concatenation is crucial for obtaining good performance in our method.

Effectiveness of h 1 (•) and h2 (•). These experiments are conducted to study the impact of h 1 (•) and h 2 (•) in our model. In the first experiment, we remove h 1 (•) from our model to evaluate the performance of h 2 (•). In this case, the output of SRU-HOS-NET before the fully connected layer is: [(vec(log( Õ1,1,2 M1 ))) T , . . . , (vec(log( Õ6,5,2 M6 ))) T ] T . In the second experiment, we remove h 2 (•) from our model to evaluate the performance of h 1 (•). The output of SRU-HOS-NET before the fully connected layer is: [(vec(log( Õ1,1,1 M1 ))) T , . . . , (vec(log( Õ6,5,1 M6 ))) T ] T . Results of these experiments are given in Tab. II. Note that the accuracy of SRU-HOS-NET degrades significantly when h 2 (•) is not used. Specifically, the accuracies of SRU-HOS-NET are decreased by 9.40, 13.09, and 17.57 percent points on DHG dataset with 14 and 28 gestures and FPHA dataset, respectively. The confusion matrices of SRU-HOS-NET on DHG dataset when using h 1 (.), h 2 (.), and both h 1 (.) and h 2 (.) are given in Fig. 5. When using only h 1 (.), SRU-HOS-NET gives more than 90% accuracy for the gestures 'grab', 'swiper', and 'swipe-+'. When using only h 2 (.), SRU-HOS-NET gives more than 90% accuracy for the gestures 'grab', 'swiper', 'swipe-+', 'rot-cw', 'rot-ccw', 'swipe v', and 'shake'. Thus, h 2 (.) captures more discriminative information than h 1 (.) for recognizing gestures with more than 90% accuracy. Moreover, when using both h 1 (.) and h 2 (.), SRU-HOS-NET gives more than 90% accuracy for the gestures 'grab', 'swipe-r', 'swipe-+', 'rot-cw', 'rot-ccw', 'swipe v', 'shake', 'expand', 'swipe-l', 'swipe-u', and 'swipe-x'. Thus, both h 1 (.) and h 2 (.) contribute to the overall performance of SRU-HOS-NET. We can also remark that the difference in performance between h 1 (.) et h 2 (.) on FPHA dataset is larger than that on DHG dataset. Note that the DHG dataset contains a set of coarse gestures that are mainly distinguished by the global movement of the hand, while all gestures of the FPHA dataset involve both the global movement of the hand and that of the fingers. This indicates that higher-order statistics could be particularly useful for finegrained gesture recognition.

C. Results on DHG dataset

The comparison of our method and state-of-the-art methods on DHG dataset is given in Tab. IV. For SRU and the method of [START_REF] Chakraborty | A Statistical Recurrent Model on the Manifold of Symmetric Positive Definite Matrices[END_REF], referred to as SRU-SPD, we ran the codes from 1,2 with the parameter settings established by the authors for obtaining their accuracies. SRU-SPD is a variant of SRU for SPD matrices, which is based on the concept of weighted Fréchet mean of SPD matrices and the Cholesky factorization for parametrization of SPD matrices. As can be observed in Tab. IV, SRU-HOS-NET and ST-TS-HGR-NET are competitive and give the best results on two evaluation protocols (14 and 28 gestures). However, in terms of network complexity, SRU-HOS-NET has far less parameters than ST-TS-HGR-NET (Tab. III). Here, the number of parameters of ST-TS-HGR-NET is about 35 times more than that of SRU-HOS-NET. Note that SRU-HOS-NET and SRU-SPD significantly outperform SRU, demonstrating the superior performance of combining first-order and second-order statistics over classical first-order statistics on this dataset. While both SRU-HOS-NET and SRU-SPD are based on the formulation of SRU, SRU-HOS-NET outperforms SRU-SPD by large margins, i.e., 8.09 and 8.69 percent points on DHG dataset with 14 and 28 gestures, respectively. This probably can be explained by the fact that (1) Our method for constructing the network input by taking into account the correlation of neighboring hand joints is effective, and (2) Similarly to SRU-SPD, SRU-HOS-NET captures first-order and second-order statistics of (c) h 1 (.) and h 2 (.)

Fig. 5: The confusion matrices of SRU-HOS-NET on DHG dataset using (a) h 1 (.), (b) h 2 (.), and (c) both h 1 (.) and h 2 (.).

Method Year Color

Depth Pose RNN/LSTM Accuracy (%) HON4D [START_REF] Oreifej | HON4D: Histogram of Oriented 4D Normals for Activity Recognition from Depth Sequences[END_REF] 2013 70.61 Novel View [START_REF] Rahmani | 3D Action Recognition from Novel Viewpoints[END_REF] 2016 69.21 1-layer LSTM [START_REF] Zhu | Cooccurrence Feature Learning for Skeleton Based Action Recognition Using Regularized Deep LSTM Networks[END_REF] 2016 78.73 2-layer LSTM [START_REF] Zhu | Cooccurrence Feature Learning for Skeleton Based Action Recognition Using Regularized Deep LSTM Networks[END_REF] 2016 80.14 Moving Pose [START_REF] Zanfir | The Moving Pose: An Efficient 3D Kinematics Descriptor for Low-Latency Action Recognition and Detection[END_REF] 2013 56.34 Lie Group [START_REF] Vemulapalli | Human Action Recognition by Representing 3D Skeletons as Points in a Lie Group[END_REF] 2014 82.69 HBRNN [START_REF] Du | Hierarchical Recurrent Neural Network for Skeleton Based Action Recognition[END_REF] 2015 77.40 Gram Matrix [START_REF] Zhang | Efficient Temporal Sequence Comparison and Classification Using Gram Matrix Embeddings on a Riemannian Manifold[END_REF] 2016 85.39 TF [START_REF] Garcia-Hernando | Transition Forests: Learning Discriminative Temporal Transitions for Action Recognition[END_REF] 2017 80.69 JOULE-color [START_REF] Hu | Jointly Learning Heterogeneous Features for RGB-D Activity Recognition[END_REF] 2015 66.78 JOULE-depth [START_REF] Hu | Jointly Learning Heterogeneous Features for RGB-D Activity Recognition[END_REF] 2015 60.17 JOULE-pose [START_REF] Hu | Jointly Learning Heterogeneous Features for RGB-D Activity Recognition[END_REF] 2015 74.60 JOULE-all [START_REF] Hu | Jointly Learning Heterogeneous Features for RGB-D Activity Recognition[END_REF] 2015 78.78 Huang et al. [START_REF] Huang | A Riemannian Network for SPD Matrix Learning[END_REF] 2017 84.35 Huang et al. [START_REF] Huang | Building Deep Networks on Grassmann Manifolds[END_REF] 2018 77.57 SRU [START_REF] Oliva | The Statistical Recurrent Unit[END_REF] 2018 72.17 SRU-SPD [START_REF] Chakraborty | A Statistical Recurrent Model on the Manifold of Symmetric Positive Definite Matrices[END_REF] 2018 78.96 ST-TS-HGR-NET [START_REF] Nguyen | A Neural Network Based on SPD Manifold Learning for Skeleton-based Hand Gesture Recognition[END_REF] hand joints' coordinates via h 1 (•). However, SRU-HOS-NET captures richer statistics than SRU-SPD thanks to h 2 (•), which can be interpreted as second-order statistics computed from h 1 (•). The method of [START_REF] Huang | A Riemannian Network for SPD Matrix Learning[END_REF] also exploits high-order statistics in the framework of CNNs. However, like SRU-SPD, it is only based on the covariance information of hand joints' coordinates. As a result, it performs significantly worse than SRU-HOS-NET, i.e., 19.16 and 19.88 percent points on DHG dataset with 14 and 28 gestures, respectively.

Fig. 1 :

 1 Fig. 1: Hand joints' positions.

9 Fig. 2 :

 92 Fig. 2: Construction of our network input.

Fig. 3 :

 3 Fig. 3: Graphical representation of our model. Solid lines indicate a dependence on the current value of a node. Dashed lines indicate a dependence on the previous value of a node. Here, T 1 = t 1 , and T 2 = t 1 + t 2 -1.

Fig. 4 :

 4 Fig. 4: The proposed network architecture. The network consists of 30 SRU-HOS, each of them produces statistics for a sub-sequence s, s = 1, . . . , 6 associated with a finger f , f = 1, . . . , 5. (a) The final outputs of each SRU-HOS are combined into Y s,f (b) Y s,f is then transformed into an Euclidean space and vectorized to give a representation of sub-sequence s associated with finger f . (c) All vectors at step (b) are concatenated to obtain a global representation of the sequence. This representation is fed to a fully connected layer and then a softmax layer. For efficient gesture recognition, we use the vectors obtained at step (c) to train a SVM classifier.

t1 y 1

 1 , 2 t1 y 2 , . . . , t1-1 t1 y t1-1 , y t1 , . . . , y M -t1+1 , t1-1 t1 y M -t1+2 , . . . , 1 t1 y M ).

  FPHA dataset. This dataset contains 1175 action videos belonging to 45 different action categories, in 3 different scenarios, and performed by 6 actors. The action categories are: charge cell phone, clean glasses, close juice bottle, close liquid soap, close milk, close peanut butter, drink mug, flip pages, flip sponge, give card, give coin, handshake, high five, light candle, open juice bottle, open letter, open liquid soap, open milk, open peanut butter, open soda can, open wallet,

  ReLU(x) = max(x, 0) is the standard linear rectifier, W r , W p , W x , W o , b r , b p , b o are the parameters of the model, R t encodes the features of averages, P t are the statistics that are dependent not only on the current input x t but also on R t , O t is the output that is fed upwards in the network, A is the set of different scales, χ χ χ t is the moving average that summarizes statistics of the sequence up to frame t. The moving averages are concatenated over all scales χ χ χ t = {χ χ χ

	(α)
	t

TABLE I :

 I Ablation analysis of our network. The best result in each row of a sub-table (t 1 , t 2 , Feature concatenation) is marked in bold.

TABLE II :

 II Effectiveness of h 1 (.) and h 2 (.). The best result in each column is marked in bold.

	Model	Number of parameters
	ST-TS-HGR-NET	672,243
	SRU-HOS-NET	18,894

TABLE III :

 III Comparison on network complexity.

TABLE IV :

 IV Performance of our method and state-of-the-art methods on DHG dataset. The best result in each column is marked in bold.

.55 14.75 11.76 1.72 2.94 1.63 1.36 72.13 5.45 1.81 3.44 8.82 4.91 1.44 4.1 1.63 83.63 1.47 1.72 1.63 86.27 1.47 1.63 2.89 2.73 8.19 3.63 85.45 1.72 1.85 1.47 1.63 5.47 1.72 3.63 7.27 87.93 3.22 3.7 1.44 3.63 93.54 12.96 1.36 3.22 75.92 1.44 1.36 3.63 1.96 3.44 77.94 4.41 88.52 1.44 1.81 5.55 78.26 5.17 1.63 1.44 100.0 10.34 1.47 1.63 11.59 84.48 1.36 1.72 82.19

	Grab	Tap	Expand	Pinch	Rot-CW	Rot-CCW	Swipe-R	Swipe-L	Swipe-U	Swipe-D	Swipe X	Swipe +	Swipe V	Shake

(a) h 1 (.)

93.1 6.55 11.76 1.72 78.68 5.45 1.96 5.88 3.27 1.63 87.27 5.88 3.44 3.27 3.63 86.27 3.63 1.85 8.19 1.36 4.91 1.81 94.54 1.72 1.61 1.63 2.73 1.72 98.27 4.83 1.85 1.47 3.27 1.44 1.36 91.93 9.25 1.44 1.81 1.61 85.18 1.81 83.82 3.27 1.72 4.91 2.94 80.32 2.89 1.85 81.15 1.72 2.89 100.0 7.24 96.55 2.89 94.52

	Grab	Tap	Expand	Pinch	Rot-CW	Rot-CCW	Swipe-R	Swipe-L	Swipe-U	Swipe-D	Swipe X	Swipe +	Swipe V	Shake

(b) h 2 (.)

94.82 8.19 1.81 13.72 3.44 78.68 1.47 6.55 92.72 2.94 1.72 3.27 1.81 86.27 1.81 8.19 1.44 1.63 96.36 1.61 1.63 1.81 100.0 1.61 1.47 3.27 1.36 1.81 96.77 3.7 1.36 94.44 4.91 1.81 92.64 1.63 3.27 1.47 78.68 95.65 1.44 100.0 100.0 1.85 1.44 97.26

TABLE V :

 V Performance of our method and state-of-the-art methods on FPHA dataset. The best result is marked in bold.

https://github.com/junieroliva/recurrent

https://github.com/zhenxingjian/SPD-SRU

D. Results on FPHA dataset

The comparison of our method and state-of-the-art methods on FPHA dataset is given in Tab. V. From this table, we can make the following observations. First, SRU-HOS-NET achieves the best result on this dataset. In particular, it outperforms the second best method ST-TS-HGR-NET by 1.39 percent point. Second, the performances of the sequential models of [START_REF] Zhu | Cooccurrence Feature Learning for Skeleton Based Action Recognition Using Regularized Deep LSTM Networks[END_REF], [START_REF] Du | Hierarchical Recurrent Neural Network for Skeleton Based Action Recognition[END_REF] are largely behind our model. The best model among them is a 2-layer LSTM, that is 14.47 percent points less accurate than our model. Note that SRU is significantly outperformed by the models of [START_REF] Zhu | Cooccurrence Feature Learning for Skeleton Based Action Recognition Using Regularized Deep LSTM Networks[END_REF], [START_REF] Du | Hierarchical Recurrent Neural Network for Skeleton Based Action Recognition[END_REF]. This is probably because those models are specially designed for skeleton-based action recognition. For example, in [START_REF] Zhu | Cooccurrence Feature Learning for Skeleton Based Action Recognition Using Regularized Deep LSTM Networks[END_REF], the authors have introduced a co-occurrence regularization term into the loss function and an in-depth dropout mechanism for the task of action recognition based on LSTMs. This result also indicates that SRU does not give good performance for hand gesture recognition on this dataset. However, the integration of high-order statistics into SRU leads to state-ofthe-art performance on this dataset. Third, compared to SRU-SPD and the networks of [START_REF] Huang | A Riemannian Network for SPD Matrix Learning[END_REF] that also exploit high-order statistics, SRU-HOS-NET significantly outperforms them by at least 10.26 percent points. Fourth, the performance gaps between SRU-HOS-NET and SRU-SPD (15.65 percent points) and between SRU-HOS-NET and SRU (22.44 percent points) are even more pronounced on this dataset.

VI. CONCLUSION

We have proposed a new RNN model for skeleton-based hand gesture recognition that integrates high-order statistics in the SRU for learning a discriminative hand gesture representation. We have evaluated the proposed method on two benchmark hand gesture datasets. On the DHG dataset, the proposed method is competitive to state-of-the-art methods, while having far less number of parameters than the best method among them. On the FPHA dataset, the proposed method outperforms state-of-the-art methods by at least 1.39 percent.