Laurent Barbiéro

Abdallahi Ould Mohamedou

Caroline Laperrousaz

Sônia Furian

Sébastien Cunnac

Polyphasic origin of salinity in the Senegal delta and middle valley

Keywords: Salinity, Soil Genesis, Aeolian Deflation, Clay Dune, Upscaling, Senegal Valley

Saline areas are a major obstacle to the development of sustainable irrigated agriculture in the Senegal valley. They have been attributed to the incorporation of marine salts in the sediments during the last marine transgression. However, this does not explain their geomorphological situation and geochemical features. They are distributed as strips about 100-200 m wide and several km long, which are composed of two parallel substrips, one located in the depressions of former creek beds, and the other on higher ground on the southern bank. In these two substrips, the chemical composition of salinity and its distribution in the soil suggest that it arises from more than one source. Comparison of saline areas of the middle valley with present-day salt accumulation in the delta suggests a 4-stage salinization process, which involves evaporation from a shallow water table, aeolian salt accumulation as clay dunes, deep transformation of shell beds into gypseous layers under temporary acid conditions induced by oxidation of pyrite under the mangrove vegetation, and secondary salinization by runoff. The 4-stage model agrees with a regional salinity chemical database and the geomorphology of the saline areas. In the Senegal valley the lower middle valley and the delta can therefore be regarded as a chronosequence, the recent salinity features occurring in the delta and the more developed ones upriver at the limit of the marine transgression.

Introduction

Salinization and alkalinization of soils and waters occur widely in arid and semi-arid regions [START_REF] Szabolcs | Salinization of soil and water and its relation to desertification[END_REF]. In the West African Sahelian zone, these problems affect all irrigation schemes, because of poor management of irrigation [START_REF] Bertrand | La dégradation des sols irrigués en Afrique de l'Ouest[END_REF]. In the Senegal river valley, located in North Senegal and South Mauritania, development of irrigated agriculture is one of the main purposes of the Manantali and Diama Dams (Organisation pour la Mise en Valeur du fleuve Sénégal [START_REF] Omvs | The objectives and the main outlines of the integrated development strategy of the Senegal river basin[END_REF]. However, the extensive saline soils, which existed prior to establishment of the irrigation schemes, are a major threat to this development. The problem is encountered especially in the delta and lower middle valley [START_REF] Loyer | Les sols salés de la basse vallée du fleuve Sénégal distribution et évolution sous cultures[END_REF][START_REF] Boivin | Soil degradation in irrigation schemes in the Senegal middle valley[END_REF].

In this region, soil salinity has been attributed to a present (in the delta) or former (in the middle valley) shallow, saline water table, resulting from a succession of recent transgressions and regressions of the Atlantic Ocean during the late Holocene (last 6000 years) [START_REF] Michel | Les bassins des fleuves Sénégal et Gambie. Etude géomorphologique[END_REF][START_REF] Brusq | Etude pédologique des cuvettes de la vallée du Lampsar (Région du fleuve Sénégal)[END_REF][START_REF] Loyer | Les sols salés de la basse vallée du fleuve Sénégal distribution et évolution sous cultures[END_REF][START_REF] Ceuppens | Impact of non drained irrigated rice cropping on soil salinization in the Senegal river delta[END_REF]. [START_REF] Michel | Les bassins des fleuves Sénégal et Gambie. Etude géomorphologique[END_REF] and [START_REF] Deckers | Major soils and their formation in the West African Sahel[END_REF] stated that the last transgression, which started about 2000 years ago [START_REF] Faure | Geoidal changes and shore-level tilt along Holocene estuaries: Senegal river area, West Africa[END_REF], was responsible for incorporation of marine salt into the sediments in which the soils have formed [START_REF] Maymard | The residual effect of submersion on soil structure (Demonstrated by applying S. Henin's index of structural instability)[END_REF]. However, preliminary studies in the middle valley showed that saline areas form roughly E-W strips >10 km long and 100-200 m wide [START_REF] Barbiéro | Salt distribution in the Senegal middle valley. Analysis of a saline structure on planned irrigation schemes from N'Galenka Creek[END_REF]. This distribution of salt does not agree with a simple scenario of salt incorporation in the sediment during the latest transgression.

In this paper, we propose a model explaining the distribution of salt, taking into account the geomorphological and geochemical features of the saline areas in the Senegal middle valley and delta, and also a regional database of salinity chemistry.

Site details

In the Senegal valley, the soils have developed from alluvial deposits following the last marine regression that occurred 2000 years ago [START_REF] Faure | Geoidal changes and shore-level tilt along Holocene estuaries: Senegal river area, West Africa[END_REF]. They form a chronosequence in which the recent soils occur in the delta and older ones, which are more strongly developed, occur upriver in the middle valley at the limit of the last marine transgression. The several processes formerly responsible for soil features in the middle valley are still active downstream in the lower middle valley or in the delta (Barbiéro et al., in press). We therefore compared salinity characteristics at two sites in the Senegal valley, one along the N'Galenka creek in the Podor region of the middle valley (northern Senegal) and the other in the Diawling National Park in the Senegal delta (southern Mauritania) (Fig. 1).

The semiarid climate of the area is characterised by a wet season (approximately 200 mm of rainfall) from July to September, a cold and dry season from October to February and a hot dry season from March to June. The temperature averages 23°C in January and 33°C in May, with daily maximum over 41°C. The prevailing winds blow from the NE (Harmattan) in the middle valley and from NE to NW in the delta because of the coastal influence. The dry and warm winds increase evaporative demand and create a plant-water imbalance. Average Class A pan evaporation exceeds average rainfall for each month and by an annual total of 2800 mm [START_REF] Verheye | Impact of climate and soil conditions on conception and implementation of irrigation schemes in the Senegal River basin[END_REF].

The middle valley

The N'Galenka creek, in the Podor region, extends from 16°26' to 16°30' N and from 14°50' to 15°05' W. Two major geomorphological units are distinguished: depressions and former river-banks (FAO-SEDAGRI, 1973), which are occupied respectively by Vertisols and Fluvisols (FAO/ISRIC/ISSS, 1998). They are known locally as Hollaldé and Fondé soils [START_REF] Maymard | The residual effect of submersion on soil structure (Demonstrated by applying S. Henin's index of structural instability)[END_REF]. Fondé soils contains about 30% clay, and their saturated hydraulic conductivity ranges from 10 -6 to 10 -5 m s -1 [START_REF] Meyer | Etude du transport des sels au cours des écoulements et de l'évolution des propriétés physiques du sol au contact de différentes solutions salines[END_REF]. Hollaldé soils contain about 60 % clay, which is composed of smectite (30 to 40 %), kaolinite (35 -40 %), interstratified illite-smectite (15 -20 %), and mica (5 -10 %) [START_REF] Mohamedou | Evolution pédogénétique des sols de mangrove de la vallée du fleuve Sénégal et phénomènes associés[END_REF]Barbiéro et al., in press). Their saturated hydraulic conductivity is around 10 -7 m s -1 . The Hollaldé soils constitute the best areas for irrigated rice growing because of their suitable physical and chemical properties. The water table occurs generally at about 7 m depth and the electrical conductivity of the groundwater ranges from 0.15 to 1 dS m -1 (Zanolin et al., 1997;[START_REF] Hammecker | Suivi et évaluation de la qualité des eaux et de l'impact environnemental des aménagements du secteur N'Gallenka[END_REF].

Surveys of soil salinity in the Podor region have been done at the paddy field scale (1-5 ha) [START_REF] De Poitevin | Etude d'impact des techniques culturales sur les aménagements hydroagricoles dans la région de Podor, Sénégal[END_REF][START_REF] Diba | Analyse comparée de deux méthodes d'évaluation de la répartition spatiale de la salinité des sols : pâte saturée et conductivimétrie électromagnétique[END_REF][START_REF] Laval | Rapport d'activités. Travaux d'expérimentation ORSTOM-FED dans la moyenne vallée du Fleuve Sénégal[END_REF] or over a larger area (288 ha) with a low density of measurements directed towards a geostatistical study of soil sampling [START_REF] Gascuel-Odoux | Variability of variograms and spatial estimates due to soil sampling : a case study[END_REF]. None of these had the objective of understanding the shape and distribution of saline areas in the landscape. More recently, an exhaustive survey of salinity covered a large area of 8000 ha at N'Galenka (Fig. 2), revealing 14 saline areas, distributed as strips that are several km long and 100-200 m wide [START_REF] Cunnac | Etude de la variabilité spatiale de la salinité des sols des futurs périmètres irrigués du Marigot N'Galenka[END_REF][START_REF] Barbiéro | Comment prévoir la répartition des sels dans la vallée du fleuve Sénégal ?[END_REF]Laperrousaz, 1999, Barbiéro et al., 2001). The strips are related to former creek channels and do not correspond to present-day salt remobilisation. A detailed study of one strip (Fig. 3) showed that it consists of two parallel substrips, a northern substrip lying in the depression of a former creek bed and a southern substrip on slightly higher ground on the southern bank of the creek [START_REF] Barbiéro | Salt distribution in the Senegal middle valley. Analysis of a saline structure on planned irrigation schemes from N'Galenka Creek[END_REF].

-4 -

The delta

In the delta, fieldwork was concentrated between the Birette and Ziré dunes close to the creek Khurumbam, in the Bell basin (Fig. 4). The main soil types here are Gypsic Fluvisols, Calcic Gypsic Fluvisols and Vertisols (FAO/ISRIC/ISSS, 1998), though Thionic Fluvisols containing pyrite (FeS 2 ) and jarosite [KFe 3 (SO 4 ) 2 (OH) 6 ] also occur near the remaining mangroves (Le [START_REF] Brusq | Nouvelles paragenèses à sulfates de fer, et de magnésium, et leur distribution dans les sols sulfatés acides du Sénégal[END_REF][START_REF] Deckers | Major soils and their formation in the West African Sahel[END_REF]. Although this site is today protected from tidal influence by a large embankment, it is partially flooded from July to September during the wet season. The water table is at 0.1 -1.2 m depth and the electrical conductivity ranges from 8 to 130 dS m -1 . Because of salinity, clay dunes are a distinctive feature of this area. From November to June, the salinity of the water and soil increases because of intense evaporation. When the topsoil dries, the increasing salinity and rapid crystallization of salts such as gypsum and halite cause the mud to flocculate forming silt and sand size aggregates, which are vulnerable to wind erosion. These aggregates accumulate on leeward sides of all wind-breaking obstacles to form saline clay dunes [START_REF] Mohamedou | Geochemistry of clay dunes and associated pan in the Senegal delta (Mauritania)[END_REF]. The environmental conditions favouring formation of clay dunes are very specific.

They require a prevailing wind during the dry season, temporary saline pan conditions, and a specific clay content in the topsoil [START_REF] Price | Physicochemical and environmental factors in clay dune genesis[END_REF]Vieillefon, 1967 ;[START_REF] Bowler | Clay dunes: their occurrence, formation and environmental significance[END_REF][START_REF] Thomas | The nature and depositional setting of arid and semi arid Kalahari sediments, southern Africa[END_REF][START_REF] Thomas | Present day lunette sediment cycling at witpan in the arid southwestern Kalahari desert[END_REF]Barbiéro et al., 1998a).

Materials and methods

Middle valley

Around Podor, six sites were selected from aerial photographs [START_REF] Geotronics | Couverture aérienne, mission Teledyne Geotronics (USA)[END_REF] in order to verify the distribution pattern of salinity. A survey of salinity was attempted by soil electromagnetic induction transects (EM38, Geonics, Ltd) [START_REF] Mcneill | Rapid accurate mapping of soil salinity by electromagnetic ground conductivity meters[END_REF] intersecting the structures suspected to be saline on the photographs. Electromagnetic conductivity (ECm) values > 150 mS m -1 were noted and saline areas were hand drawn over the background of aerial photographs. A 20 km-long saline strip, composed of two parallel substrips and crossing the whole N'Galenka region was selected for a detailed study. Four short transects across the strip (T1-T4, Fig. 2) were located using electromagnetic induction measurements to maximize the range of salinity sampled (Fig. 3). They were respectively 34, 27, 35 and 24 m long. T1, T2 and T3 cross the southern substrip of the saline area, whereas T4 crosses the northern one and is partly in the former creek bed. Soil samples along these transects were taken with a hand auger in a total of 35 profiles (8 on T1, 10 on T2 and T3, 7 on T4), each at four depths (0-0.2, 0.3-0.5, 0.7-0.9, and 1.1-1.3 m).

The electrical conductivity (EC) was measured on solutions extracted by centrifuging 1:5 soil-water mixtures using distilled-deionized water. Magnesium, sodium, calcium and potassium concentrations in the extracts were measured by atomic absorption spectrometry (Perkin Elmer).

Chloride and sulfate ion concentrations were determined by capillary electrophoresis (Waters Capillary Ion Analyser). The carbonate alkalinity was measured by HCl titration [START_REF] Gran | Determination of the equivalence point in potentiometric titrations[END_REF].

The delta

Two soil sequences were studied in the delta, across the Khurumbam creek bed and across the clay dunes fringing the creek. They were respectively 100 and 40 m long (Fig. 4). The horizons are described according to the classification for Acid Sulfate Soils of the International Institute for Land Reclamation and Improvement [START_REF] Dent | Acid sulphate soils : a baseline for research and development[END_REF]. Topography was determined using a Topochaix ® theodolite. In the first sequence, the soil morphology was reconstructed from 11 holes bored with a hand auger and 3 pits excavated on the main lateral transitions of soil horizons. In the second sequence, it was studied from 15 holes and one pit located at the top of the clay dune. 76 soil samples were taken in the second sequence reaching from in and beneath the clay dune to the creek. Saturated extracts from these samples were analysed by the methods described in 3.1. The chemical results from the second sequence were presented by [START_REF] Mohamedou | Geochemistry of clay dunes and associated pan in the Senegal delta (Mauritania)[END_REF]. Salt efflorescences were collected for SEM (XL20 Philips) observations and chemical analysis (LINK Analytical eXL energy dispersive X-ray system).

Regional data base

To estimate how representative our samples from the middle valley and delta were for developing a regional model of soil salinization, they were compared to a regional database of 1784 samples of 1:5 soil-water extracts. The data were from various locations (a-p on Fig. 1), most of them obtained before establishment of the irrigated rice growing schemes [START_REF] Maymard | Les sols du pseudo-delta du Sénégal. Mission d'aménagement du Sénégal[END_REF][START_REF] Braudeau | Etude pédologique de la cuvette de Nianga (Secteur A) -Région du fleuve Sénégal[END_REF][START_REF] Braudeau | Etude pédologique du périmètre de Balky (ex. colonat de Richard Toll), région du fleuve Sénégal[END_REF][START_REF] Loyer | Notes sur les sols et les nappes de la cuvette de N'Delle-N'Diaye[END_REF][START_REF] Loyer | Etude pédologique de la cuvette de N'Dieyrba (Région du fleuve Sénégal)[END_REF][START_REF] Brusq | Etude pédologique des cuvettes de la vallée du Lampsar (Région du fleuve Sénégal)[END_REF][START_REF] Brito | Etude pédologique du casier de Kobilo, département de Matham, région du fleuve Sénégal[END_REF]Mougenot, 1982;[START_REF] Seiny Boukar | Etude pédologique de la cuvette de N'Der, lac de Guiers[END_REF][START_REF] Loyer | Les sols salés de la basse vallée du fleuve Sénégal distribution et évolution sous cultures[END_REF][START_REF] Laval | Rapport d'activités. Travaux d'expérimentation ORSTOM-FED dans la moyenne vallée du Fleuve Sénégal[END_REF][START_REF] Samba Diene | Riziculture et dégradation des sols en vallée du fleuve Sénégal : Analyse comparée des fonctionnements hydro-salins des sols du delta et de la moyenne vallée en simple et double riziculture[END_REF][START_REF] Mohamedou | Contribution à l'étude de l'environnement du Parc National du Diawling (Mauritanie)[END_REF]Barbiéro et al., 1998b;Mohamedou et al., 2001a). Most of the analyses for this database were made at IRD laboratory (Institut de Recherche pour le Développement, previously ORSTOM) in Dakar.

Results

Salinity in the middle valley

The relationship between saline strips and geomorphological features at in several sites in the Podor region is presented in Fig. 5. The survey of soil salinity using the electromagnetic conductivity meter (EM38) reveals a unique pattern: the salt is distributed as strips along former creek beds. The saline strips are frequently divided into two parallel substrips with the northern saline substrip located in the former creek bed, and the southern fringing the creek bed on the southern bank (Figs. 5 a, c, e andf).

The chemical composition of the soil extract is plotted in Piper diagrams (Fig. 6). The composition of the least saline samples is dominated by alkalinity and Ca 2+ . In the most saline samples, the anionic composition is dominated either by Cl -or SO 4 2-, and the cations by Na + +K + .

Further discussion here is based on the distribution of Cl -and SO 4 2-

. Table 1 compares summary statistics for EC, Cl -and SO 4 2-concentrations in the soil extracts from transects at N'Galenka with the regional database. All three at N'Galenka are within the ranges for the regional database.

However, the chemistry differs because average SO 4 2-is greater than average Cl -in all four transects, but less than average Cl -in the regional database.

The concentrations of Cl -and SO 4 2-in samples from the middle valley suggest two groups of samples (Fig. 7). Group 1 is relatively scattered and the SO 4 2-concentration increases in proportion to Cl -. This group includes mainly superficial (0-0.2 m) samples and deep samples collected in the former creek bed on T4. Group 2 includes mainly deeper samples (0.3-1.3 m) dominated by SO 4 2-.

The SO 4 2-contents are fairly constant, irrespective of Cl -contents.

The distributions of EC, Cl -and SO 4 2-along the soil transects at N'Galenka are presented in Fig. 8. The contours were hand drawn using inverse distance methods between neighbouring sampling points. Maximum EC occurs around 80 cm depth in all four transects. SO 4 2-follows roughly the same distribution, but Cl -increases upwards in all the transects except at the end of T4 in the former creek bed, where maximum Cl -occurs at depth.

Soil system and chemistry in the delta

Fig. 9 shows the typical soil distribution pattern in the delta. In the creek bed at the extreme western end of sequence 1, unripe silty clay material with a high pyrite content occurs at 0.4 m depth, and is overlain by a greyish brown, silty clay, half ripe material, very acidic and with yellow mottles of jarosite. The greyish brown horizon has a coarse prismatic structure and some reddish brown (iron oxide) mottles at the base. The surface soil has a soft powdery structure.

At 0.2 m depth, near 12 m to the east, a bleached sandy horizon begins. It is discontinuous between 12 m and 24 m and its lower boundary, as seen in the pit at 24 m is abrupt, wavy and emphasized by reddish hematite mottles. The thickness of this horizon increases progressively from -7 -0.2 m to 0.95 m between 24 m and 100 m. It is locally bleached, but around the former mangrove roots it is coloured by iron oxide. Near 66 m, an abrupt change occurs in the surface structure from soft powder to centimetric platelets, which corresponds to an increase in topsoil clay content. The thickness of this clay topsoil horizon increases irregularly from 66 m to the end (100 m) of the surveyed transect, where a decimetric vertic structure appears and the clay content reaches 49 % (Mohamedou et al., 2001b).

Two layers rich in shells and shell fragments are found in the western part of sequence 1.

The first layer at 0.6 m depth is composed of Anadara senelis fragments and the second (0.9 m depth) of Thympanotonus fuscatus shells. In both layers, the shells are hard only up to 33 m and 62 m, respectively. Further east they are replaced by shells weathered by dissolution, and yet further centimetric gypseous nodules replace the two shell layers.

Aeolian accumulation of aggregates blown from the surrounding soils forms elongated clay dunes. Because the prevailing winds are from the north and north-west, the dunes are behind the remains of mangrove vegetation along the southern and south-eastern bank of the Khurumbam creek and its tributaries. They are about 1 m high, 30-40 m wide and 0.3-4 km long (Mohamedou et al., 2001b).

The soil system described in sequence 1 also occurs beneath the clay dune of sequence 2 (Fig. 9). The bleached sandy horizon is well developed but discontinuous. Two shell layers partially dissolved occur in the western part of the sequence, close to the Khurumbam. The clay dune itself is composed of stratified sandy clay loam composed of layers about 1 mm thick; 4 of them, more sandy, were identified in the pit and auger holes. Some former leaves of mangrove vegetation were observed at the base of the dune material. The floor of the clay dune is slightly undulating and higher than the surrounding land surface by about 0.4 m.

Samples taken in the delta contained greater concentrations of soluble salt than those from the middle valley. A good correlation (r 2 = 0.98, slope = 0.94) occurs between Cl -and SO 4 2-(+ in Fig. 7). In Fig. 7, the scatter plot continues that of the group 1 samples from the middle valley.

However, some points slightly departing from the correlation (asterisks in Fig. 7) are more in continuity with the scatter plot of group 2. These points are for samples taken beneath the clay dune in the shell bed and gypseous layers, and [START_REF] Mohamedou | Geochemistry of clay dunes and associated pan in the Senegal delta (Mauritania)[END_REF] have shown that the soil solution at such points is at equilibrium with respect to gypsum.

The saline efflorescences on the soil surface and in the clay dune are mainly composed of halite (NaCl) and secondarily of gypsum (CaSO 4 .2H 2 O). Mohamedou et al. (2001a) showed that their composition agrees with that of the water table. As water concentrates because of evaporation, the equilibrium with respect to these two minerals is reached, but the precipitation of thenardite (Na 2 SO 4 ) or mirabilite (Na 2 SO 4 .10H 2 O) occurs very rarely and only in the most acid pans (Mohamedou et al., 2001a).

4.3.Regional data base

The chemical results were sorted according to their chloride contents and are presented in concentration diagrams (Fig. 10). In these diagrams the concentration factor (horizontal axis) is defined as CF = [Cl]/[Cl] 0 where [ ] denotes the molarity and [Cl] 0 refers to the lowest Cl molarity measured. The concentration factor would be underestimated if halite (NaCl) precipitates, but [START_REF] Mohamedou | Geochemistry of clay dunes and associated pan in the Senegal delta (Mauritania)[END_REF] showed that even for the greatest ionic strength measured (4.43 M), the solution remains undersaturated with respect to halite. Therefore, precipitation of NaCl is not expected and chloride is a good indicator of the concentration factor. The vertical axis denotes the molarities of the major ions.

Some changes in the chemistry of the 1784 regional water samples occur along the valley.

Upriver, the soil solutions are dominated by HCO 3 -and Ca 2+ ions, as is often true in continental Africa. The more saline samples downriver are dominated by Cl -and Na + , which agrees with their marine origin. The scatter plots of Ca 2+ and SO 4 2-against Cl -concentration factor exhibit an inflection at a concentration factor between 100 and 20000, where relatively high Ca 2+ and SO 4 2occur (Figs. 10a,10b). There is a strong correlation between the CF for Cl -and Na + (r 2 = 0.99, slope = 1), although some data points depart from the straight line for CF values of 10-200 (Fig. 10c).

Discussion

Salt accumulations in the delta

Three types of salt accumulation occur at present in the Senegal delta. The first results from evaporation from the shallow water table, causing the precipitation principally of gypsum and halite in the topsoil; this will be referred to as evaporation salinity.

The second type is salt accumulation in clay dunes and will be called aeolian salinity. The clay dunes include further accumulations of gypsum and halite. The distribution of the thin sandy layers in the clay dune indicates that the dunes are fixed and grow symmetrically, which is contrary to the more common sand dunes [START_REF] Bowler | Clay dunes: their occurrence, formation and environmental significance[END_REF]. The floor of the clay dune is protected from erosion and is therefore higher than the surrounding ground, which remains exposed to wind erosion.

The third type results from formation of gypsum from biological calcite. The soil in the creek bed (west part of sequence 1 in Fig. 9) has the typical morphology of acid sulfate soil [START_REF] Dent | Acid sulphate soils : a baseline for research and development[END_REF]. At the top of the profile and eastwards along the sequence, there are standard features arising from acidity following oxidation of pyrite, and also features of acidity controlled by soil alkalinity, including the progressive dissolution of the shells and eventual disappearance of the calcareous layers. The calcium released in the saline and sulfate rich environment precipitates as centimetric nodules of gypsum, indicating the former location of the shell layers. The quantity of gypsum involved depends on the amount of shells in the sediment and the quantity of acid released by oxidation of pyrite. The third type of salt accumulation is referred to as deep gypsum salinity.

Sulfate-chloride relationships

Theory

When marine water concentrates because of evaporation, changes in the chemical composition of the solution occur by precipitation of carbonates such as calcite and sulfates such as gypsum, and possibly thenardite or mirabilite. A correlation is expected between SO 4 2-and Cl -in the concentration diagram (pathway 1-2-3 in Fig. 11). From 1 to 2, SO 4 2-increases in proportion to Cl -. Then from 2 to 3 saturation with respect to gypsum is reached and SO 4 2-increases but less rapidly than Cl -because part of the SO 4 2-is taken from the solution to form gypsum. Precipitation of sodium sulfate occurs only at greater concentrations. Inversely, if dilution occurs, the pathway is reversed and the chemical composition follows the pathway (3-2-1). From 3 to 2, the SO 4 2-content decreases less than Cl -, because some of the gypsum previously precipitated in the soil re-dissolves.

The soil solution is then at equilibrium with respect to gypsum. With further dilution (2-1), the soil solution becomes under-saturated with respect to gypsum and concentration of SO 4 2-and Cl - decrease proportionately. If dissolution of salts from aeolian salinity occurs, the relationship between Cl -and SO 4 2-will not be changed because aeolian salinity involves the same salts as evaporation salinity, i.e. gypsum and halite.

When deep gypsum occurs, the pathway followed by the soil solution during dilution in Fig. 11 coincides only partially with the theoretical evolution described above. From pathway 3 to 2, SO 4 2-decreases less than Cl -, because some of the gypsum from the evaporation salinity dissolves.

Because the solution remains at equilibrium with respect to the mineral, the deep gypsum remains undissolved. With further dilution, the soil solution becomes under-saturated and therefore aggressive with respect to gypsum, and deep gypsum begins to dissolve until saturation is reached again. The solution becomes enriched with Ca 2+ and SO 4 2-ions and evolves along pathways 4, 5 or 6 (Fig. 11); which depends on the quantity of gypsum in the soil. As long as gypsum is not exhausted, the CF for Cl -decreases but the Ca 2+ and SO 4 2-contents are maintained in the soil solution. After complete elimination of the deep gypsum, Ca 2+ and SO 4 2-decrease in proportion to the dilution effect, i.e. parallel to the slope of the lines drawn in Fig. 11.

In the samples

-10 -

The above theoretical model can explain the relationships between Cl -and SO 4 2-in the delta and in the middle valley. In the delta, the variability in the salinity is due to concentration/dilution of the soil solution, whose chemical composition possibly includes gypsum and halite from aeolian salinity dissolved during rainfall or during Senegal flooding. Dissolution of the deep gypsum is noticed only for some samples taken immediately beneath the floor of the clay dune in the gypseous layers (asterisks on Fig. 7), where the solution gets enriched in SO 4 2-and Ca 2+ . In the middle valley, the chemical composition of samples from group 1 can be related to concentration-dilution of the same soil solution, in which dissolution of salt from aeolian salinity could have occurred. The samples from group 2 with higher relative proportion of SO 4 2-and Ca 2+ cannot result only from evaporation of the same soil solution as in group 1. Nevetheless, the chemical composition in group 2 agrees with the influence of another source of salt such as deep gypsum dissolution.

Distribution of saline areas in the middle valley

The distribution of the saline strips in the middle valley is similar to that of the clay dunes in the delta. Clay dunes and saline areas are usually located along the creeks on their southern leeward bank. The saline areas have roughly similar dimensions to the clay dunes + creek beds (about 100-200 m wide and several km long), suggesting again that saline areas are related to aeolian salt accumulation in clay dunes. The strong associations between saline strips and former creek channels shown in Fig. 5 suggests that the site studied in the N'Galenka is representative of the region. Similar patterns have already been described at other sites, namely M'Boyo, Donaye and Nianga [START_REF] Barbiéro | Comment prévoir la répartition des sels dans la vallée du fleuve Sénégal ?[END_REF] around Podor. We attribute the southern substrip to the former clay dune associated with the creek and the northern substrip to the former creek bed salinized by runoff containing salt from the southern bank. The 4-stage model agrees with the morphological and geochemical features observed in the Senegal valley. Many iron pipes, which are traces of the former mangrove vegetation, are observed in the soil as far inland as Boghé, located 380 km upstream from the coast [START_REF] Michel | Les bassins des fleuves Sénégal et Gambie. Etude géomorphologique[END_REF], and attest to the previous influence of mangrove vegetation on the soil chemistry. The aeolian redistribution of salt of marine origin into clay dunes on higher ground adjacent to creek beds explains the occurrence of saline substrips on the leeward sides of previous creek channels in the middle valley [START_REF] Barbiéro | Salt distribution in the Senegal middle valley. Analysis of a saline structure on planned irrigation schemes from N'Galenka Creek[END_REF]. Both types of saline substrip have dimensions roughly similar to those of the clay dunes and creek beds. The secondary salinization of the creek bed may occur by runoff of salts from the clay dune above, as in transect T4. The model also explains the distribution of salinity chemistry in the transects, which is dominated by chloride and sulfate in the topsoil (aeolian salinity) and just sulfate at 80 cm depth (deep gypsum salinity).

Temporal model for salt evolution along the valley

The quantity of deep gypsum salinity should be influenced by variation in the initial quantities and distributions of the shells and pyrite originally present in the sediment. Isolated centimetric nodules of gypsum are common, but compact slabs of gypsum have been reported in the delta (Loyer, 1989, p. 59). In Table 1, the average composition of the transects differs from that of the regional database because desalinization is occurring in the middle valley as suggested in stage 4. At this stage, the deep gypsum is dissolving and influences the soil solution chemistry, increasing the Ca 2+ and SO 4 2-

. In contrast, in many places downriver, the soil solution is still very saline and deep gypsum dissolves very little or not at all.

The data provided by previous reports suggest the regional applicability of this scenario in the Senegal valley (Fig. 10). The straight line correlation between Cl -and Na + (Fig. 10c) from a broad range of soil solutions can be interpreted by concentration/dilution processes along the valley.

However, the organisation of the scatter plots of Ca 2+ and SO 4 2-(Fig. 10a and10b) cannot be explained only by concentration/dilution, but with the model including the dissolution of deep gypsum salinity. Our conclusions agree with those of Samba Diene (1998), who attributed the accumulation of deep gypsum to a process independent of the present soil solution chemistry at a site in the lower middle valley.

On Fig 10c, the samples that do not fit the strong correlation between sodium and chloride molarities for concentration factors between 10 and 300 suggest that small quantities of sodium sulfate (thenardite or mirabilite) are possibly incorporated into the deep gypsum, and dissolve with it when dilution occurs. Such precipitation could occur under very acidic conditions, as shown by Mohamedou et al. (2001a). These samples were collected in a former clay dune at M'Boyo, located about 20 km upstream from Podor. The unusual chemistry of the saline soil solution may be specific to this location or may predict a change in the chemical composition of the saline areas upriver from Podor. We do not have pertinent data to explain the deviation of these samples from the correlation more definitely.

Conclusions

Previous researches in the middle valley of the Senegal River have attributed the formation of saline areas to marine salts deposited in the sediment during the most recent marine transgression. However, this fails to explain the present distribution and features of the saline areas.

Based on their morphology and geochemistry, the presence of saline areas results from the presence of former clay dunes, but their chemical characteristics are explained by a 4-stage salt accumulation involving: 1) concentration in the groundwater because of evaporation, 2) aeolian deflation and clay dune formation, 3) transformation of shell banks into gypseous layers under acidic conditions generated by oxidation of pyrite, and 4) local salinization by runoff. The 4-stage accumulation of salts has developed under restricted environmental conditions, which are similar to those observed today in the delta, where their evolution continues. The regional applicability of this model is supported by a database collected along the valley.

Our results suggest that the soils in the valley form a chronosequence, in which the most recent occur in the delta and the earlier, more developed, soils occur upriver near the limit of the last marine transgression. We have already used the chronosequence concept to describe the pedological processes involved in the development of soil profiles in the valley (Barbiero et al., in press). It may also be useful in studying the evolution of other soil features, such as the natural control of the acidity of the mangrove swamps and the influence of water alkalinity on possible soil degradation resulting from irrigation. -31 - Table 1. Comparison of the descriptive statistics of the EC (mS cm -1 ), Cl and SO 4 contents (meq l -1 ) measured on the soil extracts in transects T1-T4 in the N'Galenka area with the regional database of 1784 samples.

  These observations suggest a new 4-stage model explaining the distribution of salinity and the variable composition of the soil solution (Fig.12). Stage 1 was the period of accumulation of pyrite by sulphur-fixing bacteria in a waterlogged environment supplied with seawater sulfate under a tidal influence. In Stage 2, the drainage resulting from the marine regression brought oxygen into the previously waterlogged soils, and the pyrite was oxidized to form sulfuric acid. At this stage, only the creeks were still occupied by mangrove vegetation. The shell beds dissolved and became gypseous horizons. The shallow water table evaporated and the rapid crystallization of gypsum and halite in the topsoil caused the development of a soft powdery structure. The aggregates of the powdery structure were vulnerable to wind deflation, and therefore accumulated behind wind-breaking obstacles as clay dunes, especially on the southern banks of the creeks. The first two stages are currently observed in the Senegal delta.In Stage 3, the water table was at a lower level and no longer subject to evaporation. The rains of the wet season and the seasonal overflows of the Senegal River resulted in partial desalinization of the soils. Evaporation salinity was eliminated from areas around the clay dunes by the percolating water, but the clay dunes themselves remained saline for a longer time because of the greater quantity of salt in them. By dissolving some superficial gypsum, percolating water became saturated with respect to this mineral and therefore non-aggressive on reaching the deep gypsum. The deep gypsum began to dissolve only after almost complete elimination of the superficial gypsum (Stage 4). Surface runoff could then cause secondary salinization in depressions located near former clay dunes, i.e. the creek bed. Where secondary salinization occurred, the chemical composition of the salinity in the creek bed would have been similar to that of the superficial part of the soil. Conversely the composition of the topsoil salinity differed from the deep gypseous horizons, because chloride content was greater in the superficial layer.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Location of the studied sites in the Senegal valley. a-p are sites where the geochemical database was collected. a: Khurumbam; b: Keur Macene; c: N'Diaye; d: Gorom; e: Lampsar; f: Boundoum; g: N'Der; h: Rosso; i: Balky; j: Dagana; k: Fanaye; l: Nianga; m: N'Gaoule; n: Donaye; o: M'Boyo; p: Kobilo.

Fig. 3 :

 3 Fig. 3: Topography (a) and distribution of electromagnetic conductivity (b) in a saline strip showingthat the northern substrip is located in a depression and the southern substrip on higher ground (modified fromBarbiero et al., 2001). (.) denotes measurement points.

Fig. 5 :

 5 Fig. 5: Aerial photographs showing strip pattern of saline areas and their relationship to former creek channels at six sites in the Podor region.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Piper diagrams showing the variability in the chemical composition of the soil aqueous extracts along the transects studied in the N'Galenka area. Open circles denote the 5 least saline samples.

Fig. 8 :

 8 Fig. 8: Distribution of electrical conductivity, chloride and sulfate contents along vertical cross sections of transects T1-T4 in the N'Galenka area ( ▪ denotes sampling points).

Fig. 9 :

 9 Fig. 9: Soil sequences 1 and 2 across the Khurumbam creek and clay dune in the Diawling park ( ▪ denotes sampling points on sequence 2).

Fig. 10 :

 10 Fig. 10: Relationships between sulfate (a), calcium (b), sodium (c) and a concentration factor derived from chloride contents along the Senegal valley. Letters in (c) and the respective concentration ranges refer to the locations in Fig. 1.

Fig. 12 :

 12 Fig. 12: The 4-stage model for the origin of saline areas in the Senegal valley.

Acknowledgement

This work was partially realised in the framework of the PSI-CORAF (Regional Coordination of the Irrigated Systems Pole of the Conference of Agronomic Research Directors in Africa). We thank Samba Seck and Alfred Savadogo for help in collecting data, and Dr. Vasanthi Dass for editorial advice.

The salinization process is currently occurring only in the delta and no longer in the middle valley, where only the salt present at about 60 cm depth can be remobilised by field irrigation. As currently observed in many irrigation schemes, this can cause problems in agriculture. Because the presence of saline areas cannot be easily detected from the soil surface, rapid mapping of saline soil distribution using electromagnetic induction methods and aerial photographs to identify former geomorphological features, is an essential prerequisite for the establishment of new irrigation schemes in this region.


Antoine Blais

Nicolas Couellan
email: nicolas.couellan@recherche.enac.fr

Evgenii Munin

A novel image representation of GNSS correlation for deep learning multipath detection

Keywords: Deep learning, GNSS, Multipath, Convolutional neural networks, Correlation

This paper proposes a novel framework for multipath prediction in Global Navigation Satellite System (GNSS) signals. The method extends from dataset generation to deep learning inference through Convolutional Neural Network (CNN). The process starts at the output of the correlation stage of the GNSS receiver. Correlations of the received signal with a local replica over a (Doppler shift, propagation delay)-grid are mapped into grey scale 2D images. They depict the received information possibly contaminated by multipath propagation. The images feed a CNN for automatic feature construction and multipath pattern detection. The issue of unavailability of a large amount of supervised data required for CNN training has been overcome by the development of a synthetic data generator. It implements a well-established and documented theoretical model. A comparison of synthetic data with real samples is proposed. The complete framework is tested for various signal characteristics and algorithm parameters. The prediction accuracy does not fall below 93% for C/N0 ratio as low as 36 dBHz, corresponding to poor receiving conditions. In addition, the model turns out to be robust to the reduction of image resolution. Its performance is also measured and compared with an alternative Support Vector Machines (SVM) technique. The results show the undeniable superiority of the proposed CNN algorithm over the SVM benchmark.

Introduction

The dissemination of GNSS receivers in smartphones and cars has made it natural to each of us to have access to our localisation at any time. However, the quality of the position calculated by a GNSS equipment may be reduced when the received signal is degraded. This degradation can nd its origin in a defect of the signal generation system, carried by the satellite, it is the evil waveform case. The receiving conditions can also be a source of disturbance, this is typically the case when interferences or multipaths are in addition to the useful signal.

More specically, inside the GNSS chip, the calculation process of the Position Time Velocity (PVT) solution relies on the accurate synchronisation of the receiver on the signal transmitted by each satellite in view. This is achieved by correlating the received signal with a local replica signal which parameters are under control. The correlation is done against the three unknown parameters of the incoming signal which are required to achieve synchronisation: the propagation delay τ , the Doppler shift δf (due to the relative motion between the satellite and the receiver) and the carrier phase θ. The aim of the receiver is to nd the set of estimators which maximizes the correlation, resulting in the tightest synchronisation of the local replica signal. Any alteration to the signal of interest will skew the estimators and biased the position delivered to the user.

A large amount of research and analysis has been conducted so far to detect, classify, identify and nally mitigate these impairments. As the GNSS receiver has to track the direct signal by mean of a Delay-Locked Loop (DLL) to estimate the propagation delay, multiple methods have been proposed which use the already existing correlator outputs required by this DLL. The narrow correlator technique [START_REF] Van Dierendonck | Theory and performance of narrow correlator spacing in a GPS receiver[END_REF], the early-late-slope technique [START_REF] Townsend | A practical approach to the reduction of pseudorange multipath errors in a Ll GPS receiver[END_REF], the strobe correlator [START_REF] Garin | Enhanced strobe correlator multipath rejection for code & carrier[END_REF], the double-delta correlator [START_REF] Mcgraw | GNSS multipath mitigation using gated and high resolution correlator concepts[END_REF] and the multipath intensive delay lock loop [START_REF] Jardak | Multipath insensitive delay lock loop in GNSS receivers[END_REF] are among the most representative methods of this class. They all take advantage of the geometric shape of the auto-correlation function of the Pseudo-Random Noise (PRN) code, as dened later on in (2) and (3) and illustrated in Figure 3, to detect and mitigate the multipath distortion. Their relative simplicity is their principal benet at the expense of their eectiveness. On the other hand, more sophisticated techniques, yet demanding in hardware resources, have been developed. In the statistical approach, the Multipath Estimating Delay Lock Loop (MEDLL) is a reference implementation of the maximum likelihood principle [START_REF] Townsend | Performance evaluation of the multipath estimating delay lock loop[END_REF]. It matches the correlator outputs with candidates of multipath auto-correlation functions parameterized by magnitudes, delays, and phases. The shortest estimated delay is then retained as the one of the direct path. The frequency domain has also been explored, through the Fourier transform (Zhang & Bartone, 2004a) or the wavelet decomposition [START_REF] Zhang | Real-time multipath mitigation with wavesmoothtm technique using wavelets[END_REF]. Indeed, due to their specic spectral characteristics in comparison to the direct path, the multipath can be identied and excised. However, these methods may damage the signal of interest, especially when the multipath frequencies are close to the spectrum of the direct path.

To overcome the limitations of these classical signal processing methods, Machine Learning (ML) techniques have also been considered. Starting from the early 2000s, some research work has been dedicated to the use of ML techniques to facilitate the error mitigation in GNSS signals. For instance, an hybrid neural network architecture based on multilayer perceptron to mitigate multipath error for Low Earth Orbit (LEO) satellites has been proposed [START_REF] Vigneau | Neural networks algorithms prototyping to mitigate GNSS multipath for LEO positioning applications[END_REF]. Later, with the advances of kernel methods, the authors of [START_REF] Phan | A unied framework for GPS code and carrier-phase multipath mitigation using support vector regression[END_REF] were able to develop a support vector regressor to mitigate multipath on ground xed Global Positioning System (GPS) stations and using signal geometrical features. Other similar studies were conducted with various choices of features construction. For example, in [START_REF] Hsu | GNSS multipath detection using a machine learning approach[END_REF] and [START_REF] Suzuki | NLOS multipath classication of GNSS signal correlation output using machine learning[END_REF], nonline of sight (NLOS) multipath detection is carried out using features directly extracted from the correlator output.

The recent and signicant advances in Articial Intelligence (AI), and notably in ML, have opened up new perspectives. In [START_REF] Quan | Convolutional neural network based multipath detection method for static and kinematic GPS high precision positioning[END_REF], using a CNN, a carrier-phase multipath detection model is developed. The authors propose to extract feature map from multi-variable time series at the output of the signal processing stage using 1-dimensional convolutional layers. Deep learning spoong attack detection in GNSS systems was addressed in the research literature [START_REF] Schmidt | A GPS spoong detection and classication correlator-based technique using the LASSO[END_REF] as well. Handcrafted features based on early-late phase, delay and signal level from the correlation output of the tracking loop were used to train a deep fully-connected neural model. A review of the recent applications of ML in GNSS is also proposed in [START_REF] Siemuri | Machine learning utilization in gnssuse cases, challenges and future applications[END_REF], focusing on use cases relevant to the GNSS community.

The method proposed in this article aims at making use of an ecient CNN architecture for multipath detection. The intent is to exploit the full power of CNN by letting the convolutional mechanism construct its own feature space from the whole correlator information. Indeed, features are not extracted from the signal but the signal is rather transformed into 2D images in the timefrequency domain. No correlator output signal information is lost during the process and the CNN is able to build its own representation of corrupted/non corrupted correlated signals.

The search ranges for the values of the propagation delay and the Doppler shift spans a 2D grid which forms in turn 2D-images at the output of the correlation process, the correlator output in short. Regarding the phase value, the phase estimation error is captured over [0, 2π] by mean of two orthogonal projections. These projections are carried out by the correlation with the In-phase (I) signal replica on one side and the in-Quadrature (Q) signal on the other. A diagram representing this process is given in Figure 1. As depicted, the correlation operation is implemented through a product followed by an integrate and dump stage and generates 2D image representations of the I and Q channels. These images coded into 3D tensors will feed a downstream CNN, as it will be seen below.

This study proposes a complete framework to train and assess a CNN model on correlator output 2D-images in order to detect whether the GNSS signal is subject to multipath or not. A graphical representation of this framework is given in Figure 2. Our technique exploits the full power of deep learning architectures by sampling the complete correlated signal information in the timefrequency domain and the I and Q channels. Features are not handcrafted but rather constructed automatically by the convolution mechanism that elab-

Receiver antenna port X X I Q r ( t ) s (t ; τ ,δ f , θ ) s (t ; τ ,δ f , θ + π / 2) Integrate & Dump ∫ T i Integrate & Dump ∫ T i
Local replica signal generator orates its own representation of the relevant feature space to detect multipath corrupted signals.

The main contributions of this study can be listed as follows:

-Raw and complete information from the GNSS correlator outputs are synthesized in 2D-images. The correlation delay and Doppler shift ranges are selected in order to capture complete multipath information. This is a novel approach in comparison to standard multipath mitigation techniques that are using only one dimensional delay correlation information.

-A CNN model is used to automatically extract relevant features for multipath detection from the images of correlator outputs. This contributes to the very recent and emergent use of modern ML techniques in the GNSS signal processing eld.

-The proposed framework covers the generation of image data, the choice of the CNN architecture, its training as well as its validation. Experiments are fully reproducible. To the best of our knowledge such complete workbench is unique in the research community.

-The achieved average detection accuracy for realistic multipath parameters ranges in standard receiving conditions is above 93%. This performance has been shown to be robust to the reduction of correlator output image resolution.

The organization of this article follows the framework depicted in Figure 2. In Section 2, the GNSS signal model used in this work is presented, then the correlation process is introduced and a model of its output is detailed, for both the direct path and the multipath signal. Section 3 exposes the dataset elaboration using a specic software generator. Next, in Section 4, the CNN model proposed in this paper to detect the multipath contamination is explained in detail. In Section 5, the experiments conducted to assess the proposed technique are exposed and the results discussed. Finally, Section 6 draws conclusions and perspectives from this work.

Problem statement 2.1. GNSS signal model

The fundamental principle behind the calculation of the user position by a GNSS receiver is trilateration. It implies the measurement of the geometric distances between the antenna of the receiver and satellites of known positions. Indeed, a distance d in particular is estimated through the propagation delay aecting the signal during its propagation from the satellite to the receiver antenna, τ = d/c with c the speed of light. This is made possible by a specic signal structure, recalled in equation (1) which models the signal at the antenna port [START_REF] Won | Signal processing[END_REF]:

r(t) = √ 2CD(t -τ )c(t -τ ) cos(2π(f c + δf )t + θ) + b(t) (1)
where

• C is the power of the received signal,

• D(t) is the navigation message, binary encoded (±1),

• c(t) is the PRN code sequence, specic to each satellite,

• f c is the carrier frequency,

• b(t) is an Additive White Gaussian Noise (AWGN) which accounts for the thermal noise of the receiver, referred to the antenna port.

The results presented in this paper where established using the PRN code sequences of the GPS L1 C/A legacy signal. However, the authors are condent that they could be generalized to other navigation signals, with the same structure, as no specic assumption has been made on c(t).

In this model, the receiving condition of a signal in particular is assessed by its C/N0 gure, in other words the ratio of the signal power to the Power Spectral Density (PSD) level N0 of the (white) noise b(t). Clearly, the accuracy of the estimation of the related distance d will depend upon this C/N0 ratio. Without any other perturbation than the noise, the quality of the nal position calculated by the receiver, from a set of distances d at its disposal, is then completely determined by the corresponding set of C/N0 ratios, along with the relative satellites-receiver geometry though.

A model of the I and Q correlator outputs [START_REF] Won | Signal processing[END_REF]) is as follows:

I = AK(∆τ ) cos(π∆f T i + ∆θ) sinc(π∆f T i ) + n I (2) Q = -AK(∆τ ) sin(π∆f T i + ∆θ) sinc(π∆f T i ) + n Q (3) with • T i the integration time,
• A a coecient depending on C, D and T i ,

• ∆τ the propagation delay estimation error,

• ∆f the Doppler shift estimation error,

• ∆θ the phase estimation error,

• K(∆τ ) the auto-correlation function of the PRN code in ∆τ ,

• n I and n Q the noise components.

It is worth noting the sinc function behavior of the correlator outputs as a function of the Doppler shift estimation error.

Figure 3 gives a graphical representation of the noise-free I and Q correlator outputs as functions of ∆τ and ∆f , for ∆θ = 0. The shape of the autocorrelation function of the PRN code can be observed in the ∆f = 0 plane, for delays in units of µs. The previously mentioned sinc function is visible in the ∆τ = 0 plane, for Doppler shifts in units of Hz.

Multipath contamination

As some perturbations can distort the desired signal, the received signal can not always be modeled simply using equation (1). Among these perturbations, multipath is considered to be an important source of degradation [START_REF] Kaplan | Understanding GPS : principles and applications[END_REF]. This is especially the case in urban environment, inducing reduced positioning accuracy. Multipath is due to the reection of the direct signal path on a surface in view of the receiver. As a consequence, a specic multipath can be modeled in the same way as the direct signal in (1):

m(t) = 2C MP D(t -τ MP )c(t -τ MP ) cos(2π(f c + δf MP )t + θ MP ) (4)
where C MP , τ MP , δf MP and θ MP have the same denition as in Section 2.1, but for the multipath. Due to the larger propagation distance of the multipath in particular, it is to be noted that C MP ≤ C and τ MP > τ . What is more, depending on the time-varying relative geometry of the satellite-receiver-reector system, there is no reason for δf MP being equal to δf nor θ MP having the same value as θ.

In general, a receiver is impacted by multiple multipaths, especially in urban environments where reectors are numerous. Sometimes, the direct path may even be absent due to an obstruction, for example when high buildings are surrounding the receiver [START_REF] Ziedan | Improved multipath and NLOS signals identication in urban environments[END_REF]. However, in this study the assumption is made that the direct path is always present and a single multipath will be considered.

Synthetic dataset generation

Availability of GNSS multipath data

In order to test our prediction models, an articial signal generator was developed. The data are generated in the form of two matrices, one for each of the I and Q channels, according to equations ( 2) and (3). The axes of these matrices are in Doppler shift estimation error ∆f and code delay estimation error ∆τ . The output data corresponding to this main signal can be parameterized as a function of the coherent integration time T i in ms and the carrier-to-noise ratio C/N0 in dBHz.

Noise sample production

At the correlator output the noise is not only spatially correlated inside each I and Q image, but also cross-correlated between them. The exact derivation of the auto-correlation and cross-correlation functions of the noise are still to Figure 4: Empirical auto-correlation and cross-correlation functions of the noise at the correlator output for the I and Q channels, PRN number set to 10 establish. To overcome this impossibility to generate the noise contribution at the correlator output from an analytical model, a workaround has been developed. A signal r(t) (1) made of a simple noise term b(t) is correlated according to the process described in Figure 1, as would be a true received signal. This correlation process is implemented in a software GNSS receiver developed by the SIGnal processing and NAVigation (SIGNAV) research team of the ENAC laboratory. The noise samples available at the correlator output are then collected and stored in a dataset, to be added on demand as n I and n Q in equations ( 2) and (3). Figure 4 gives an example of empirical auto-correlation and crosscorrelation functions of noise samples for PRN number 10.

An illustration of the noisy output of the synthetic data generator for T i = 20 ms and C/N0 = 45 dBHz is given in Figure 5, with the corresponding attened images in Figure 6 (note the value of the navigation bit D = -1 this time, in comparison to Figure 3 where D = +1).

Assessment of the synthetic data

In order to validate the synthetic data generator its outputs have been compared to the data from an IFEN SX3 GNSS receiver. Two dierent data collection sessions have been conducted.

1. The receiver has been fed with a signal produced by a Spirent GSS6560 generator. The scenario implemented in the generator simulates the takeo and initial climb of a commercial aircraft from runway 14L of the Toulouse-Blagnac Airport. The ight happens on Tuesday the 28th of May 2019 from 12:55 UTC. The multipaths are disabled in this scenario. 2. A high end GNSS antenna has been connected to the receiver. The antenna was set up in a clear view site to avoid multipath contamination. In both cases the SX3 receiver sampling frequency was set to 20 MHz. The samples were stored for post-processing by the software GNSS receiver already mentioned in section 3.2. The resulting reference images are available on [START_REF] Blais | A synthetic GNSS correlator output generator[END_REF] for the two sessions.

Figures 7 and8 provide an example of visual comparison of images. It is worth noting that these real data were used only to validate the generator. Indeed, for training purposes it would be unrealistic to produce a sucient amount of labeled physical signals.

Model for GNSS multipath data generation

The considered model integrates I and Q signals, I and Q multipath and correlated receiver noise. If a multipath signal is received in addition to the main signal, as the signal processing chain is linear, the correlator output can then be considered as the sum of the correlator output of the main signal and the one due to the multipath. In this work, a single multipath contamination is considered. Its contribution (I MP , Q MP ) to the correlator output is considered as an additional term to the main signal (I, Q), detailed in (2) and (3):

I = I + I MP (α MP , ∆τ MP , ∆f MP , ∆θ MP )
(5)

Q = Q + Q MP (α MP , ∆τ MP , ∆f MP , ∆θ MP ) (6)
where

• α MP = C MP /C < 1
is the multipath attenuation coecient in comparison to the main path, • ∆τ MP = τ MP -τ > 0 is the code delay in excess to the main signal delay,

• ∆f MP = δf MP -δf is the dierence between the Doppler shift of the main signal and the multipath,

• ∆θ MP = θ MP -θ is the dierence between the phase of the main signal and the multipath.

The I/Q image dataset generation

The data generator has been implemented with the Python language [START_REF] Blais | A synthetic GNSS correlator output generator[END_REF] to produce datasets of I and Q images according to the signal and noise models detailed in Sections 3.2 and 3.4. This software is referred as the generator in this paper. The generator is fully congurable with respect to the following parameters:

• α MP , ∆τ MP , ∆f MP and ∆θ MP which entirely dene the multipath. The denition intervals of these parameters are detailed in Appendix A. Their probability distributions are claried in Section 5 dedicated to the experiments conducted with the help of the generator,

• C/N0 ratio which sets the strength of the direct path signal with respect to the receiver noise,

• N the number of pixels along each of the delay and Doppler shift axes. The denition domains of the delay and Doppler shift axes are set in Appendix B. The size of the images is then N × N pixels. It is important to mention that 2N 2 is then the number of correlators required to implement the technique proposed in this paper. Hence, it is a direct measure of its complexity as the correlation operation is from far the most power and time consuming process in a GNSS receiver,

• T i the coherent integration time. In this study T i = 20 ms, a value corresponding to the duration of one navigation bit D as dened in equation ( 1). It ensures the longest correlation time, so the best accuracy for τ , δf and θ estimation, without bit transition during the correlation which would otherwise lower the nal result.

The experiments led in this work were all done with datasets of 600 pairs of I and Q images, equally split in 300 with multipath and 300 without. However, the generator can provide datasets of arbitrary size and distribution on demand.

Convolutional neural network model 4.1. Image classication using CNN

Convolutional neural networks [START_REF] Goodfellow | Deep Learning[END_REF] are nowadays considered as among the most powerful tools to learn information from images. This is explained by their computational eciency and their impressive performance on image information processing. Their learning capability comes from their ability to automatically construct and combine abstract features from an image. The rst layers of the network are composed of convolutional layers. Such layers apply several lters over the various regions of the input image and create feature maps that are various versions of the ltered input image. Several convolutional layers are usually stacked in order to progressively extract meaningful information from the feature maps as the depth of the network is increased. The last layers of the CNN perform the task of classication usually through several layers of neurons with dense connectivity. The weights of each convolution layer lters (also called kernels) and those of the dense layers are learned through supervised learning based on gradient back-propagation. The underlying structure of convolutional layers has the advantage of having sparse connectivity and high weight sharing among neurons, which leads to much greater computational efciency over fully connected neural networks when image sizes are large or the image resolutions are high.

The above principles are at the heart of CNN architectures. Several additional components are usually integrated in the network. Some Pooling layers are used in order to reduce the dimension of feature maps. A Flatten layer is added in front of the stack of dense layers so as to transform feature map information in vector-like input. To increase the generalization power of the network, a dropout mechanism can also be used to articially and randomly remove a small portion of the neuron connection within the network. The number and the organization of these various components generate several possible CNN architectures [START_REF] Shin | Deep convolutional neural networks for computer-aided detection: CNN architectures, dataset characteristics and transfer learning[END_REF].

Choice of a CNN architecture

Among CNN architectures that have proven to be eective in practice [START_REF] Voulodimos | Deep learning for computer vision: A brief review[END_REF], the VGG-like architecture has been shown to be one of the best choice for image feature extraction [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF]. Despite Figure 9: The CNN architecture used in this study: the input has 2 channels corresponding to I and Q channels, the rst and second convolutional layers have 16 and 32 lters respectively with ReLu activations and the pooling layer is a 2 × 2 max pooling operation layer (this gure has been generated by the NN-SVG tool (LeNail, 2019)). its number of parameters to be trained compared to other popular and more complex architectures such as Inception V3 [START_REF] Szegedy | Rethinking the inception architecture for computer vision[END_REF], ResNet [START_REF] He | Deep residual learning for image recognition[END_REF] and other variants, it has been widely adopted in practice by the ML community. The architecture is composed of several blocks of convolutional layers that are each separated by a pooling layer that decreases the feature map dimension between blocks. As the depth increases, the dimension of layer input decreases but the number of lters increases. When data are organized on a multi-scale basis, mixing macro and micro patterns, the number of convolutional blocks should be chosen suciently large. For the GNSS multipath application, such multi-scale representation is not expected in the signal. This is the reason why the chosen architecture only includes one convolutional block composed of two convolutional layers. It is therefore a very simple instance of a VGG-like network.

VGG architectures are usually implemented for RGB images (meaning three input channels: the 'R', 'G' and 'B' channels). For the specic case of I/Q images, each input image is actually composed of two channels (I channel + Q channel) that are sharing the GNSS correlation signal information. Input images are therefore tensors of size N × N and depth 2. Figure 9 provides the precise network architecture and layer dimension used in this study. 4.3. Feature maps and heatmaps as a combined mean to visualize the feature space Although articial neural networks are usually thought of as blackbox models that construct complex decision surfaces that are dicult to interpret for a specic application, their CNN instances provide some ways to get knowledge on the features that are designed during the training process. The convolutional Figure 10: Examples of feature maps extracted from the last convolutional layer for the generated synthetic data using 3 × 3 lters. layers extract meaningful information in the image for the learning task at hand. Therefore, visualizing the constructed lters after each convolutional layers gives information on the structure of the image that was extracted by the layer, the so-called feature maps. Visualizing the stack of the various feature maps at the output of the various layers will provide insightful knowledge for the user, see Figure 10 for an example of such stack of lters constructed by the proposed CNN. The gure shows ltered GNSS correlation images from the last convolutional layers. It provides understanding on the regions of the signal image that are extracted for multipath detection. Clearly the various peaks of the signal are considered for the detection task.

In addition, there is also the possibility to build a class-discriminative localization map that will highlight the regions of the input image that are most important for the trained CNN to assign a specic class to the image. Several methods have been proposed along this line [START_REF] Zhou | Learning deep features for discriminative localization[END_REF][START_REF] Selvaraju | Grad-CAM: Visual explanations from deep networks via gradientbased localization[END_REF]. The GradCAM method [START_REF] Selvaraju | Grad-CAM: Visual explanations from deep networks via gradientbased localization[END_REF], that will be used later in this study, is often used for such purpose as it is applicable to most standard CNN architectures. The main idea is to compute the gradients of the score y c of a specic class c with respect to the feature maps A k ∈ R u×v of a convolutional layer. A class importance weight α c k is then constructed by averaging these gradients over the height and width of the feature map as follows:

α c k = 1 u × v i j ∂y c ∂A k ij .
These weights are then used to build a weighted combination of activation maps during a forward propagation of a specic input image:

F c GradCAM = ReLU k α c k A k .
Using the ReLU function will account only for positive inuence of class c. To visualize the corresponding inuence in the input space, resizing will be necessary as the activation map F c usually does not have the same dimension as the input. In the experiments, an example of such GradCam activation map is provided, see Figure 12.

In Section 5.7, such heatmaps are constructed for the multipath detection task in the I/Q images. These visualization tools provide further understanding on the patterns of the correlation signals that are important for multipath discrimination. This information may be used later to further rene the signal frame that carries relevant information.

Experiments

Experimental setup

This section describes the experiments which were conducted to evaluate the performance of the proposed model on the dataset described in Section 3.5. For each test case the mean, median and standard deviation values of accuracy and F1 score averaged over 20 runs are provided. The range and probability distribution of the multipath parameters as dened in equations ( 5) and ( 6) are set, unless otherwise specied, as follows:

• α MP the multipath attenuation coecient is uniformly distributed in [0.1, 0.9],

• ∆τ MP the additional propagation delay of the multipath is also uniformly distributed, in [0, 3T c /2] as explained in Appendix A,

• ∆f MP the dierence in Doppler shift between the direct signal and the multipath is distributed in [-125, +125] Hz according to a truncated zeromean normal distribution with a standard deviation set to 125/3,

• ∆θ MP the dierence between the phase of the main signal and the multipath is uniformly distributed in [0, 2π].

The value of the navigation bit D, as dened in equation ( 1), is chosen randomly with equal probability between -1 and +1 for each pair of I and Q images.

Inuence of multipath characteristics

In this section, the inuence of the multipath characteristics on the detection performance is assessed. The experimental conditions are set in this way: C/N0 = 47 dBHz and the image resolution is 80x80. This corresponds to good observation conditions for the direct path, so that the CNN response can be clearly observed. The tests were conducted on 5 equally spaced values for both the dierence in Doppler shift ∆f MP from 0 to 50 Hz and the additional propagation delay of the multipath ∆τ MP from 0 to T c . The results are represented in Tables 1 and2. They show that there is no signicant inuence of the Doppler shift on the detector performance. On the other hand, when the propagation delay approaches 0.2T c (meaning that the multipath gets close to the main signal), a slight degradation of the accuracy and F1-score by 6% can be observed.

Inuence of the signal to noise ratio C/N0

In this section the inuence of C/N 0 ratio on the detection performance is assessed. The experimental conditions are the following: T i = 20 ms and the image resolution is 80x80. The multipath parameters distributions are identical to those dened in Section 5.2. The tests were conducted on equally spaced values of C/N 0 ratio from 24 to 46 dBHz. The results are presented on Table 3. From this results the current model shows high robustness towards the noisy incoming images for values down to 36-38 dBHz. Then, as expected, the model performance decreases greatly with the value of the C/N 0 (when I and Q images are noisier). On Figure 11, it can be also observed that the standard deviation of the metrics decreases when the image becomes less noisy. This means that the detection model achieves greater robustness in performance when the noise decreases. A value of around C/N 0 = 36 dBHz seems to be also a threshold after which the performance is much higher.

Inuence of the I/Q image resolution

In order to evaluate the performance of the multipath detector, the proposed algorithm was applied on images of various resolutions N ∈ {20, 40, 60, 80} with the intent to estimate the best compromise performance vs image resolution. The multipath parameters distributions are the same as before. The results in Table 4 show that the model performance decreases with the image resolution. However, for resolutions above 40 pixels per axes, the model classication performance stays above 95%.

Discussion on CNN performance results

The experiments exposed in this section put in light that the performance of the proposed method is more than adequate in the three axis of importance • The multipath parameters have been varied over their respective realistic range to assess the validation accuracy. It is not less than 94% in any case,

• The C/N0 ratio, that reports the receiving condition of the GNSS signal of interest, was changed from 46 (good) to 24 (poor) dBHz. A failover is observed in the performance around 36 (fairly poor) dBHz, which establishes a quite acceptable operational limit to our network,

• The image resolution, parametrized by N the number of points per axes, has been gradually decreased to measure the robustness of the algorithm to hardware limitations. Indeed, a validation accuracy of 93% is still achieved for N = 20.

With respect to the chosen architecture, the experiments conrm that the detection task at hand does not require large depth as not much multi-scale learning seems to be needed. Indeed, our choice of a quite shallow architecture has proven to be appropriate for ecient multipath detection.

Comparison with an SVM based detection technique

In this experiment, the performance of the proposed method is compared to the SVM multipath detection technique proposed in [START_REF] Suzuki | NLOS multipath detection by using machine learning in urban environments[END_REF]. The SVM model constructs a maximum margin separating hyperplane between the two classes of data points (multipath/no multipath). SVM are often used as an alternative method to neural network as they are able to handle non linearly separable data (Vapnik, 1998;[START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF]. Since the method proposed in [START_REF] Suzuki | NLOS multipath detection by using machine learning in urban environments[END_REF][START_REF] Suzuki | NLOS multipath classication of GNSS signal correlation output using machine learning[END_REF]) also collects signal information from the output of the correlator block, from a data collection point of view, it is similar to the technique proposed in this study. However, unlike CNN, they require prior data feature engineering. For comparison purposes, the feature extraction pipeline used in [START_REF] Suzuki | NLOS multipath detection by using machine learning in urban environments[END_REF] was also implemented on 13 correlator outputs and identical SVM hyperparameters were used. More specically, the features proposed in [START_REF] Suzuki | NLOS multipath detection by using machine learning in urban environments[END_REF] were extracted as follows:

• Number of local maxima of the correlation outputs per period

F 2 = N local-maxima ∆t
where ∆t is the correlation interval taken equal to coherent integration period.

• Distribution of the delay of the maximum correlation output However, in [START_REF] Suzuki | NLOS multipath detection by using machine learning in urban environments[END_REF], the authors have also used a signal strength versus elevation angle feature (referred as F 1 in their article). This feature was not taken into account here as no physical context was introduced in the experiments (receiver's speed, satellite constellation) since the generated dataset is synthetic. In the experiment, the C/N 0 is taken high (47 dBHz) for both algorithms to assess performance in favorable noise conditions. Table 5 reports the results on the mean, median and standard deviation values of accuracy and F1 score1 averaged over 20 runs of both models (SVM and CNN). The results show that the proposed CNN algorithm considerably outperforms the SVM benchmark.

F 3 = 1 M M i=1 (t i-max -t) 2

Analysis of CNN automatic feature construction: towards CNN multipath detection interpretability

The experiments have highlighted the multipath detection power of the CNN when compared to SVM. The automatic feature extraction that takes place in the convolutional layers is able to catch the geometrical dependencies in the data. To further demonstrate this property, activation maps as described in Section 4.3 have been computed using the GradCam technique (see Figure 12). Clearly, the activation maps show hot detection regions around both signal and multipath peaks (red and yellow areas on the gure) while cold regions that do not carry detection relevant information are marked in dark blue. Distortion around the main signal is highlighted in the activation map in various multipath situations (near or far from main signal peak). When there is no multipath, no distortion appears in the heatmap and a much more symmetrical pattern can be observed. These observations in a sense validates the detection mechanism and provides a sound and clear interpretation of the CNN decision rule.

Conclusions

In this study, a complete GNSS multipath detection framework based on deep learning has been presented. The proposed method starts with the construction of training image data from synthetic receiver correlation outputs. A precise image generation process based on specic parameter denition intervals is described. This process optimizes the relevance of built-in information within dataset samples. A CNN architecture is then presented and tested with the constructed dataset. For various multipath parameter choices, experiments have demonstrated the detection performance of the proposed deep learning model. Further investigation using heatmaps provides additional understanding of the detection model decision rule and validates its relevance. The results provided in this study are very encouraging and should motivate further research combining ML techniques and GNSS signal processing modelling. More specically future research will focus on multiple multipath that characterizes urban environment. Investigations will also be conducted with deep regression architectures for multipath parameters estimation. The time dynamic of the multipath should also be studied in order to improve current static learning models. For the Doppler shift observation interval, Figure A.14 shows that the upper bound can be taken |δ f | < 5.5/T i because for higher frequencies the value of the multipath sinc function becomes negligible in comparison to its maximum. Nevertheless, as in the previous section, depending on the value of T i , |δ f | is

  Synoptic view of the correlation process. The received signal is correlated with two local replica signals in quadrature whose parameters span a grid. The two correlator outputs form then 2D-images which fed a downstream CNN. The tilde notation indicates the local parameter by opposition to the received signal unknown parameter.

Figure 2 :

 2 Figure 2: Proposed framework to train and assess a CNN model on correlator output 2Dimages (with references to article sections).

Figure 3 :

 3 Figure 3: An illustration of the noise-free I and Q correlator output model, ∆θ set to 0, PRN number to 1

Figure 5 :

 5 Figure 5: An illustration of the I and Q correlator outputs of the synthetic data generator, PRN number set to 1

Figure 8 :

 8 Figure 8: Comparison between an output of the synthetic data generator and a real sample from a SX3 receiver, Q channel, PRN number set to 1

Table 3 :

 3 Figure 11: Average validation accuracy (left) and average F1 score (right) with respect to carrier to noise ratio C/N 0 (vertical bars represent standard deviation values)

Figure 12 :

 12 Figure 12: Example of heatmaps (bottom) with respect to input image sample (top)

Figure A. 14 :

 14 Figure A.14: Limit case of the dierence in Doppler shift between the direct path and the multipath

  Figure A.13 represents the correlation between the local signal and the sum of the direct path and the multipath signals as a function of the delay oset between the local and received signals. The triangle pattern corresponds to K(∆τ ) (2,3), the specic auto-correlation function of the PRN code, with T c the bit period of the code. This corresponds to the limit case when the multipath signal arrives with the same Doppler shift and with the same amplitude and phase as the direct path signal. The two circles at -T c /2 and +T c /2 represent the two extreme correlation points used in a classical GNSS receiver to track the signal of interest. The case without disturbance of these measurement points is realized when ∆τ MP > T c + T c /2, so ∆τ MP is taken in [0, 3T c /2].Figure A.14 represents the correlation between the local signal and the sum of the direct path and the multipath signals as a function of the frequency oset between the local and received signals. This corresponds to the hypothetical limit case when the multipath signal arrives with the same code delay and with the same amplitude and phase as the direct path signal. The disturbance induced in the correlation function of the direct path by the multipath sinc function is considered to be negligible for |∆f MP | > 2.5/T i as it peaks under 10% of the value at ∆f = 0. However, depending on the value of T i , |∆f MP | is further bounded by the physical value of 800 Hz. This is the maximum Doppler shift which could be encountered for example by a commercial aircraft in the approach or landing phase at a speed of 140 knots. Finally, |∆f MP | < min(800, 2.5/T i )

Table 1 :

 1 Mean (µ) and standard deviation (σ) of prediction accuracy and F1-score with respect

	to Doppler shift
	∆f dop (Hz)	0	10	20	30	40	50
	µ Acc 0.97 0.96 0.96 0.98 0.97 0.99
	σ Acc 0.01 0.01 0.02 0.02 0.01 0.01
	µ F 1 0.96 0.99 0.97 0.97 0.97 0.99
	σ F 1 0.01 0.02 0.02 0.03 0.02 0.01
	Table 2: Mean (µ) and standard deviation (σ) of prediction accuracy and F1-score with respect to propagation delay
	∆τ (T c )	0.2	0.4	0.6		0.8 1.0
	µ Acc	0.94	0.94	0.98		0.96 1.0
	σ Acc	0.01 0.023 0.022 0.012 0.0
	µ F 1	0.94	0.95	0.98		0.95 1.0
	σ F 1 0.008	0.02 0.023 0.011 0.0

Table 5

 5 where t i-max is the code delay of the maximum correlation output, t is the mean of the code delay, and M is the number of correlator output samples.

	: Comparison of Mean (µ) and standard deviation (σ) of prediction accuracy and F1-score for SVM and CNN models Model SVM CNN
	µ Acc	0.73	0.98
	µ F 1	0.74	0.95

The F1 score is not to be confused with the F 1 feature from(Suzuki et al.,

2017).

Acknowledgements

This project has been partly funded by the SESAR Joint Undertaking under the European Union's Horizon 2020 research and innovation programme under grant agreement No 783287. The opinions expressed herein reect the authors' view only. Under no circumstances shall the SESAR Joint Undertaking be responsible for any use that may be made of the information contained herein.

The authors are greatful to the Direction Technique de l'Innovation (DTI) of the French Civil Aviation Authority (DGAC) for providing a simulation environment to validate our synthetic generator.

Appendix A. Denition intervals of the multipath parameters

In this section the ranges of variation of the multipath parameters ∆τ MP and ∆f MP are addressed.