Cécile Bessou
email: cecile.bessou@cirad.fr

Laurence D C Chase

Ian E Henson

Amir F N Abdul-Manan

Llorenç Milà

I Canals

Fahmuddin Agus

Mukesh Sharma

Melissa Chin

Cécile Bessou

Pilot application of PalmGHG, the RSPO greenhouse gas calculator for oil palm products

Keywords: Palm Oil, GHG, Calculator, RSPO, Scenario testing

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Palm oil is now the most used vegetable oil worldwide, accounting in 2011 for 31.3% of the global oils and fats production (Oil World, 2012). About 10% of global production is certified by RSPO, the Roundtable on Sustainable Palm Oil (Oil World, 2012;[START_REF] Rspo | RSPO manual on best management practices (BMP) for existing oil palm[END_REF]. RSPO is a non-profit association created in 2003 and legally registered in 2004. It promotes the production and consumption of sustainable palm oil through a voluntary certification scheme. This certification scheme was developed by the stakeholders of the seven sectors involved in the palm oil commodity chain: growers, processors/traders, social NGOs, environmental NGOs, manufacturers, retailers, and banks. It relies on the growers and processors for compliance with 39 principles and criteria 1 (P&Cs) of sustainability that were defined by consensus in 2007 2 . RSPO is directed by an executive board (EB), the highest RSPO authority, which plays a major role in organising an annual general assembly and supervising member groups. During 2009-2011, the EB commissioned two successive science-based working groups on greenhouse gas (GHG WG) tasked with the aim of identifying ways leading to meaningful and verifiable reductions in GHG emissions.

Within the framework of P&Cs revision, the GHG WG2 recommended to the EB the use of PalmGHG [START_REF] Chase | PalmGHG Calculator: A greenhouse gas accounting tool for palm products, Beta-version. The Roundtable for Sustainable Palm Oil -RSPO[END_REF], a GHG calculator that allows producers to calculate the GHG balances of oil palm products. PalmGHG was developed by the GHG WG2 as an Excel spreadsheet using the life cycle assessment approach. It is based on a tool developed previously by [START_REF] Chase | A detailed greenhouse gas budget for palm oil production[END_REF]. PalmGHG quantifies the major sources of emissions and sequestration for a palm oil mill and its supply base, and is thus consistent with the operation of the certification scheme. The calculator is flexible, allowing for different crop rotation lengths and use of alternative default values. It calculates the total net emissions per ha, allocates these to products, and expresses them as t CO 2 e/t palm product, e.g. crude palm oil (CPO). The calculations are done annually: allowing for identification of principal emission sources for management purposes and for regular reporting, and communication across the supply chain and with interested parties. This paper presents briefly the structure and content of PalmGHG Beta version 1.a (of December 2012), and focuses on the results from the pilot phase of PalmGHG carried out in 2011 with nine RSPO companies.

Identified needs for further research and PalmGHG improvements are also stressed.

Material and methods

1 The 39 RSPO P&Cs cover 8 dimensions:1) Commitment to transparency [2 criteria], 2) Compliance with applicable laws and regulations [3 criteria], 3) Commitment to long-term economic and financial viability [1 criterion], 4) Use of appropriate best practices by growers and millers [8 criteria], 5) Environmental responsibility and conservation of natural resources and biodiversity [6 criteria], 6) Responsible consideration for employees and for individuals and communities affected by growers and mills [11 criteria], 7) Responsible development of new plantings [7 criteria], and 8) Commitment to continuous improvement [1 criterion]. 2 In 2012 the first revision of the P&C took place (this process is planned for every 5 years); at the moment this manuscript was written the proposed revision of the P&C was due to be voted in April 2013. PalmGHG is a significant development of the GWAPP3 model of [START_REF] Chase | A detailed greenhouse gas budget for palm oil production[END_REF] which, while straightforward to use, has a number of limitations in that palms are assumed to have an even age distribution, no allowance is made for variations in inputs and outputs within the palm area, and several aspects involve modelling rather than being based on field data. These limitations were addressed and the model simplified by assessing only the main sources of emissions and sequestration, and by using as much field data as possible, resulting in a smaller, more flexible, and more site-specific calculator. The following description of PalmGHG content is relevant for PalmGHG Beta version 1.a (of December 20124 ).

PalmGHG provides an estimate of the net GHG emissions produced during the palm oil production chain.

The emissions are presented as CO 2 equivalents (CO 2 e). Following the IPCC guidelines (IPCC, 2006), the GHGs considered are CO 2 , N 2 O, and CH 4 . Conversion factors to give CO 2 e (global warming potentials) are those determined by IPCC (2007), and correspond to a 100 year timeframe; i.e. 298 kg CO 2 e/kg N 2 O and 25 kg CO 2 e/kg CH 4 . The conversion factor for biogenic CH 4 is calculated from the ratio of the molecular weights of CO 2 and CH 4 to account for the released CO 2 originating from photosynthetic fixation; i.e. a global warming potential of 22.25 kg CO 2 e/kg CH 4 [START_REF] Wicke | Different palm oil production systems for energy purposes and their greenhouse gas implications[END_REF].

Net GHG emissions are calculated by adding the emissions released during land clearing, crop production and crop processing, and subtracting from these emissions the sequestration of carbon in the standing crop and conservation areas. The system boundary is shown in Figure 1. In the first step, net emissions are calculated as tonnes of CO 2 e per hectare. From the yield of Fresh Fruit Bunches (FFB) and the extraction rates in the mill, product-specific emissions are calculated per tonne of Crude Palm Oil (CPO) and per tonne of Palm Kernel (PK). Allocation of the net emissions of CO 2 e between the various products of interest CPO and PK, and subsequently between Palm Kernel Oil (PKO) and Palm Kernel Expeller (PKE), is carried out according to the extraction rates of each mill. However, allocation factors based on the energy contents of oils and coproducts may be also provided in future programmed versions of PalmGHG, in order notably to compile PalmGHG results with further assessments up to biodiesel production within the RED (European Commission, 2009). System expansion is otherwise parameterised for valued co-products getting out of the system, i.e. excess electricity from burning fibre and shell or methane capture. No emission burdens are calculated for co-products recycled within the system. However, should these internally recycled co-products be exported, further system expansion scenarios should be implemented.

The emission sources included in the calculator are: i) Land clearing;

ii) Manufacturing, transport and use of fertilisers;

iii) Nitrous oxide and carbon dioxide resulting from the field application of fertilisers and mill coproducts; iv) Fossil fuel used in the field, mainly for harvesting and collection of FFB; v) Fossil fuel used at the mill; vi) Methane produced from palm oil mill effluent (POME); and vii) Carbon dioxide and nitrous oxide resulting from the cultivation of peat soils.

In addition, the following GHG fixation and carbon credits are also considered: i) Carbon dioxide fixed by oil palm stand, ground cover and in plantation litter;

ii) Carbon dioxide fixed by biomass in conservation areas (methodology still under development and discussion); iii) GHG avoided by the selling of mill energy co-products (e.g. electricity sold to the grid; palm kernel shell sold to industrial furnaces).

These ten items account for the bulk of the GHG emissions and fixation occurring during the oil palm crop cycle [START_REF] Chase | A detailed greenhouse gas budget for palm oil production[END_REF]. Items that are not included in the budget are the nursery stage, pesticide use, fuel used for land clearing, emissions embodied in infrastructures and machines, and the sequestration of carbon in palm products and co-products. These items are generally negligible GHG sources or sinks [START_REF] Schmidt | Life cycle inventory of rapeseed oil and palm oil[END_REF][START_REF] Choo | Determination of GHG contributions by subsystems in the oil palm supply chain using the LCA approach[END_REF]. The other emissions are particularly small when annualised over the crop cycle.

Provision is made for separate budgets for a mill's own crop (produced on estates that belong to the same company as the mill) and outgrower crop (such as produced by independent or contracted farmers who are smallholders). PalmGHG uses the annualised emission and sequestration data to estimate the net GHG balance for the palm products from both own and outgrower crops at an individual mill. Emissions from the bi-omass cleared at the beginning of the crop cycle are averaged over the cycle. Emissions from the other sources are averaged over the three years up to and including the reporting date, thus simplifying data collection and smoothing out short-term annual fluctuations. The estimates can be updated on a yearly basis to reflect changes in operating conditions and the growth and age distribution of palms.

INSERT FIGURE 1 ABOUT HERE

Scientific background

The scientific background is detailed in [START_REF] Chase | PalmGHG Calculator: A greenhouse gas accounting tool for palm products, Beta-version. The Roundtable for Sustainable Palm Oil -RSPO[END_REF], but some insights on the main assumptions are given in this paper. Emphasis was put on using the most updated equations and data relevant for palm oil production systems and tropical agricultural conditions.

Land use change

The approach used to evaluate the contribution of land use and land use change (LULUC) to net GHG emissions is based on the Stock-Difference Method of IPCC (2006) Tier 1. The time interval to assess the stock-difference corresponds to the average crop cycle length. The latter is defined by the user and can differ between own-crop and out-grower areas, as well as between plantations on mineral or peat soils. It usually slightly varies around 25 years, which is set as default value. Carbon stocks considered in the various land uses include both above-and below-ground biomass. Changes in soil organic matter in mineral soils might be significant in the long term but were not considered as the evidence concerning such change is limited and contradictory especially in the long term [START_REF] Soussana | Carbon cycling and sequestration opportunities in temperate grasslands[END_REF][START_REF] Seguin | Moderating the impact of agriculture on climate[END_REF][START_REF] Germon | Évaluation des impacts environnementaux des Techniques Culturales Sans Labour en France[END_REF].

Areas planted each year with a breakdown of previous land uses are recorded. Based on literature data, values for eight previous land uses and oil palm stands are currently available in PalmGHG (Table 1). Additional previous land uses, such as primary forest or shrub, should be implemented soon following a review process by experts under the supervision of RSPO GHG secretariat. User-defined previous land uses may also be implemented provided evidence of field measurements and a validation by experts based on a dedicated guidance by RSPO GHG secretariat. Emissions arising from land clearing are calculated based on measured carbon contents or in their absence on an assumed carbon content of 45% in the biomass of the previous vegetation.

In agreement with Stock-Difference Method of IPCC ( 2006), carbon stock in crop stand is accounted for, i.e.

the incremental carbon that is fixed in the oil palm stand over the whole crop cycle of about 25 years. It corresponds to some threshold also provided in PAS2050 (BSI, 2011). Carbon fixed in palm products (oils) and co-products (fibre, shell etc.) is too short-lived to be accounted for (BSI, 2011). Given the emphasis on gathering accurate data on oil palm agronomy, carbon sequestration in the oil palm stand are modelled with OPRODSIM 5 and OPCABSIM 6 models [START_REF] Henson | OPRODSIM, a versatile, mechanistic simulation model of oil palm dry matter production and yield[END_REF][START_REF] Henson | Modelling carbon sequestration and greenhouse gas emissions associated with oil palm cultivation and land-use change in Malaysia. A re-evaluation and a computer model[END_REF], which are specifically designed to estimate oil palm stand biomass largely based on Malaysian conditions. The total amount of carbon sequestered in the reporting year is calculated by multiplying the area of each year of planting by the amount of carbon sequestered, adding these together, and dividing by the total area to give t C/ha/yr. Field observations revealed that biomass growth and yields are generally lower in the case of out-growers [START_REF] Chase | A detailed greenhouse gas budget for palm oil production[END_REF][START_REF] Khasanah | The carbon footprint of Indonesian palm oil production[END_REF]. To reflect this difference, contrasting simulation scenarios of crop sequestration can be used as default estimates within PalmGHG for mill own crops and out-growers. A "vigorous growth" simulation model is considered for own crops, and an "average growth" simulation is used for out-growers. Alternative models exist [START_REF] Van Noordwijk | Estimating the carbon footprint of biofuel production from oil palm: Methodology and results from two pilot areas in Indonesia[END_REF]updated by Khasanah et al. 2012;[START_REF] Asmara | Oil palm plantation carbon stock calculator[END_REF]. These models, as well as potential allometric measurements carried out by the user, may be further implemented in future versions of PalmGHG to account better for management impacts on the biomass stand.

Field emissions

Emissions due to fertilisers contribute significantly to total agricultural GHG emissions and so affect the final GHG balance of palm oil [START_REF] Yusoff | Feasibility Study of Performing a Life Cycle Assessment on Crude Palm Oil Production in Malaysia (9 pp)[END_REF][START_REF] Pleanjai | Greenhouse gas emissions from the production and use of palm methyl ester in Thailand[END_REF][START_REF] Arvidsson | Life cycle assessment of hydrotreated vegetable oil from rape, oil palm and Jatropha[END_REF][START_REF] Choo | Determination of GHG contributions by subsystems in the oil palm supply chain using the LCA approach[END_REF]. Therefore, they have been accorded special attention in PalmGHG. N 2 O direct and indirect field emissions are calculated according to IPCC Tier 1 (IPCC, 2006). Following the same guidelines (IPCC, 2006), CO 2 emissions from urea application which is subject to volatilisation losses are also accounted for.

Provision is given for nine widely used synthetic fertilisers and two organic ones (Empty Fruit Bunches, EFB, and POME) but additional fertiliser types can be included by the users if required. Methane emissions due to POME are accounted for at the mill stage and vary according to the type of treatment. Calculations of CH 4 production and amounts and losses during digestion, flaring, or electricity production are based on factors from [START_REF] Schmidt | Life cycle inventory of rapeseed oil and palm oil[END_REF] and the Environment Agency (2002). Radiative forcing of CH 4 from POME accounts for its biogenic origin [START_REF] Wicke | Different palm oil production systems for energy purposes and their greenhouse gas implications[END_REF] and is set to zero when CH 4 is flared and converted to CO 2 .

Finally, both CO 2 and N 2 O emissions related to peat soil cultivation are included in the assessment. N 2 O emission factor is taken from Tier 1 IPCC (2006), i.e. 16 kg N-N 2 O/ha.year. CO 2 emissions are based on the results of the RSPO Peatland Working Group, which placed emphasis on the importance of managing the water table depth to limit peat soil subsidence and related CO 2 emissions. CO 2 emissions due to peat cultivation are hence currently calculated using Eq. 1 (Hooijer et al., 2010), and vary with water table management. However this may be revised in the future to better account for further management practices and better differentiate peat oxidation and root respiration. In PalmGHG a default drainage depth when water table is actively managed is set to 60 cm, considered, as good management practice (RSPO, 2012) to be the maximum 5 Oil Palm Production Simulator 6 Oil Palm Carbon Budget Simulator 7 water level to be maintained below the peat surface. From Eq. 1 this leads to annual emissions of 54.6 t CO 2 /ha. Peat CO 2 emission (t CO 2 /ha.year) = 0.91 x Drainage depth (cm)

Eq. 1

Further research is needed to better define how agricultural management and in particular water table management might influence the amount of CO 2 and N 2 O emissions linked to peat cultivation [START_REF] Rspo | RSPO manual on best management practices (BMP) for existing oil palm[END_REF].

PalmGHG is a flexible tool and options are given to change default parameters. For instance, the crop cycle length can be adapted. It can notably be changed to reflect differences between crops on mineral soils and those on peat soil, which are often shorter due to accentuated sensitivity to pest and diseases and poor anchorage of palm roots (Wetlands International, 2010).

INSERT TABLE 1 ABOUT HERE

Mill emissions

At the mill level, two main sources of GHG emissions are present, fossil fuel consumption and methane emission from POME, although only the latter is significant at the supply chain level [START_REF] Pleanjai | Greenhouse gas emissions from the production and use of palm methyl ester in Thailand[END_REF][START_REF] Choo | Determination of GHG contributions by subsystems in the oil palm supply chain using the LCA approach[END_REF]. POME emissions vary depending on its treatment, i.e. conventional digestion, flaring or conversion to electricity. Parameters to account for the impact of the various treatment options are provided based on [START_REF] Yacob | Baseline study of methane emission from open digesting tanks of palm oil mill effluent treatment[END_REF][START_REF] Yacob | Baseline study of methane emission from anaerobic ponds of palm oil mill effluent treatment[END_REF][START_REF] Schmidt | Life cycle inventory of rapeseed oil and palm oil[END_REF] Given the field variability of some of these key parameters, e.g. the amount of POME per tonne of fresh fruit bunches processed or the amount of methane emitted per tonne of POME treated, further development of PalmGHG might include new parameters to better define POME emissions according to detailed treatment practices (including potential co-composting).

The pilot process

PalmGHG is a GHG calculation tool that aims to help palm oil producers to identify GHG hotspots as well as management options to reduce these emissions. To this end, the content and format of the tool were discussed by scientists and supply chain stakeholders for about two years. We reached compromises combining consistent scientific background with practical knowledge from the producers. The system boundary encompasses the whole supply area of one mill, which is the certification unit and hence the management unit too.

The format of input data was also agreed in order to ease data collection and yearly update, while seeking for a sound smoothing-out of intrinsic variability. Default values are provided in order to ensure the consistency of the calculation, but options are given to the user to implement actual data, which shall help to better account for actual management options. As much as possible, parameters linking detailed practices with emission determinism were introduced in order to identify management control levers. Finally, the results are ex-pressed in disaggregated units in order to manage GHG emissions at several levels, while avoiding problem shifting.

A pilot study was carried out with nine RSPO companies. This pilot phase was crucial both to check for the consistency of the tool, especially in terms of data availability, and to determine the ease of use and suitability of PalmGHG as a management tool. In June 2011, a preliminary questionnaire was sent to correspondents from the pilot companies. This questionnaire was the starting point of correspondence between these companies and RSPO appointees who were responsible for guiding company correspondents in the use of PalmGHG. The questionnaire aimed to introduce PalmGHG to the correspondents and help them to define the boundaries for the chosen mills and identify the data needs. Mail exchanges, as well as field visits, assisted the compilation of input data and calculation of GHG balances. The correspondents and RSPO appointees worked in parallel on PalmGHG with the same data sets to allow for cross-checking and trouble shooting.

Results were presented to all the stakeholders during the annual RSPO conference in November 2011.

Results

PalmGHG pilot results

Results from ten mills are presented in this paper (Table 2). The nine companies involved in the pilot provided data for several of their mills so that complete data sets were finally consolidated for ten mills in total.

The studied mills cover the main worldwide producing countries, i.e. Indonesia, Malaysia, as well as Brazil

and Papua New Guinea. The average GHG balance was 1.67 t CO 2 e/t CPO, and ranged from -0.02 to +8.32t CO 2 e/t CPO. Main hotspots of emissions were land clearing, peat oxidation, fertiliser-related emissions (embodied emissions in fertiliser production and field emissions) and methane from POME. Across the mills without supply from peat area, land clearing represented 41-80% of total emissions, POME represented 15-35% of total emissions, and fertiliser-related emissions represented 3-19% of total emissions. Field emissions are a common hotspot in GHG assessment of agricultural productions due to fertiliser use, and adding land clearing emissions reinforces the overall weight of the field emissions. In the case of palm oil, these two hotspots are even more important that field operations are mostly done by hand (e.g. organic fertiliser application, harvest). In the pilot mills, fossil fuel used in the field only contributed to 0-5% of total emissions.

Fuel use at the mill stage also was very limited (0-2% of total emissions) due to the recycling of co-products to produce heat and power. The main GHG source at the mill stage was methane emissions from POME. At the time of the assessment, none of the pilot mills was capturing methane from the POME digestion to flare it or convert it to electricity. Currently, only 10% of the mills in Malaysia and Indonesia have the installation to capture methane from effluent digestion (EPA, 2012). The contribution of the various emission sources and sink are shown for mill H in Figure 2.

In the cases where fruits were supplied from peat land area, CO 2 and N 2 O emissions due to peat cultivation became overwhelming. For the two mills encompassing significant peat area in their supply area, CO 2 and N 2 O from peat cultivation accounted for 62-77% and 8-10% of total emissions; respectively. Extreme values of the observed range of GHG balances are hence explained by peat emissions on the one end, and by land use change on the other end. Indeed, negative GHG balances can be reached if biomass stand in previous land uses was lower than that of palm oil stand, while other sources of emissions remain relatively low (i.e.

notably without peat emissions). Among the pilot mills, which all had identical POME treatment, previous land use and the percentage of the area on peat soil were the main causes of GHG balance variation.

INSERT TABLE 2 ABOUT HERE INSERT FIGURE 2 ABOUT HERE

Data for outgrowers were available for five mills. On a hectare basis, contribution of outgrowers to annual supply varied between 2-29% with a median value of 7%. Fruit yields are generally lower in plantations of outgrowers than own crops due to non-selected planting material and less efficient fertiliser management [START_REF] Bessou | LCA of Palm Oil in Sumatra, Comparison of Cropping Systems[END_REF][START_REF] Harsono | Energy balances and greenhouse gas emissions of palm oil biodiesel in Indonesia[END_REF]. In the pilot mills, yields of outgrowers' plantations were between 4-60% lower than those of own crops. Impact on the GHG balance at the field level was however more affected by the differences in previous land uses and the proportion of peat land in own crops (Fig. 3). In the pilot, outgrowers were not cultivating peat soils and primarily planted palms after arable crops or grassland.

INSERT FIGURE 3 The contribution of mill fuel is negligible and not visible on the graph. Net emissions below 0.3 t CO2e/t CPO can thus be obtained from a mature industry that is replanting palms, capturing methane and generating electricity from the biogas (Scenario 5 in Figure 4). This was highlighted in recommendations of the RSPO GHG WG2 to the RSPO EB within the frame of P&Cs revisions.

INSERT FIGURE 4 ABOUT HERE

Feedback from users

The main difficulties encountered by the pilot companies were related to data collection. Some difficulties may diminish as the recording process becomes routine. However, collecting data from outgrowers may remain a problem unless they are guided in the processes of data collection and PalmGHG use.

Correspondents from the pilot companies were reactive and provided feedback on the design and the ease of use of PalmGHG. Some modifications have already been implemented, but the following points still need further consideration: 1) use of graphs to help visualise the effects of adjusting key variables; 2) allowing for a more immediate way to change crop lifespans; 3) incorporation of further routines to calculate the standing biomass of crops and natural vegetation; 4) incorporation of additional fertiliser types such as compounds and mixtures; and 5) identification of further reliable peat emission parameters. These further developments will be made possible through both software programming and continuous research effort on improving the science basis behind emission factors and process modelling.

Discussion

GHG balances calculated with PalmGHG are within the range of those found in the literature. However, depending on the system boundaries and particularly on assumptions regarding land clearing and peat emissions, estimated GHG balances vary greatly around 2.3 t CO 2 e/t CPO [START_REF] Schmidt | Life cycle inventory of rapeseed oil and palm oil[END_REF], from 0.6-1 t CO 2 e/t CPO [START_REF] Siangjaeo | Implications of land use change on the life cycle greenhouse gas emissions from palm biodiesel production in Thailand[END_REF], or 2.8-19.8 t CO 2 e/t CPO [START_REF] Reijnders | Palm oil and the emission of carbon-based greenhouse gases[END_REF]. In previous studies, the relative importance of different sources of GHG emission was shown to be similar; land clearing and the cultivation of peat soils are the two most important contributors [START_REF] Germer | Estimation of the impact of oil palm plantation establishment on greenhouse gas balance[END_REF][START_REF] Reijnders | Palm oil and the emission of carbon-based greenhouse gases[END_REF][START_REF] Wicke | Different palm oil production systems for energy purposes and their greenhouse gas implications[END_REF]. Some studies that do not directly address them still mention the primary importance of these factors [START_REF] Yusoff | Feasibility Study of Performing a Life Cycle Assessment on Crude Palm Oil Production in Malaysia (9 pp)[END_REF][START_REF] Pleanjai | Greenhouse gas emissions from the production and use of palm methyl ester in Thailand[END_REF][START_REF] Stichnothe | Life cycle assessment of two palm oil production systems[END_REF]. In all studies methane from POME and fertiliser production and use are also important contributors [START_REF] Choo | Determination of GHG contributions by subsystems in the oil palm supply chain using the LCA approach[END_REF][START_REF] Pleanjai | Greenhouse gas emissions from the production and use of palm methyl ester in Thailand[END_REF][START_REF] Siangjaeo | Implications of land use change on the life cycle greenhouse gas emissions from palm biodiesel production in Thailand[END_REF], although their relative importance depends on whether land use change and peat emissions are included . As shown in the PalmGHG scenario tests, it is often emphasised that methane capture can result in significant GHG reductions in comparison with open-pond POME treatment, e.g. of between 30 to 50% [START_REF] Vijaya | Life cycle inventory of the production of crude palm oil -A gate to gate case study of 12 palm oil mills[END_REF][START_REF] Chuchuoy | Preliminary Study of Thailand Carbon Reduction Label: A Case Study of Crude Palm Oil Production[END_REF].

Carbon stocks in biomass stands and peat emissions are very sensitive parameters. Research efforts are still needed to better quantify carbon stocks and impacts of agricultural practices on these stocks, especially in the case of peat cultivation. As highlighted by [START_REF] Couwenberg | Greenhouse gas emissions from managed peat soils: is the IPCC reporting guidance realistic?[END_REF], IPCC coefficients for GHG emissions from managed peat soils may be still too low. Thanks to the work of the RSPO Peatland Working Group, some management factor could be introduced within PalmGHG to account for the impact of the depth of drainage on peat oxidation. It already allows for the producers to see the impact of peat management on GHG emissions in order to foster best practices implementation. However, further improvement is needed to better model the impact of peat management on the various GHG emissions. This is of paramount importance in Southeast Asia where peat lands represents 57% of the tropical peat area and 77% of the tropical peat carbon pool (10-14% of the global peat carbon pool). A great part of tropical peat land area is furthermore located in the two top palm oil producing countries, Indonesia with 21 million ha (65%) and Malaysia with 2.6 million ha (10%) [START_REF] Page | Global and regional importance of the tropical peatland carbon pool[END_REF].

We acknowledge that the modelling of LULUC impact on the GHG balance is not sensitive to the actual residence time of carbon in the various land uses. The debate around timing of emissions and removals in LCA has been on-going for the last twenty years or so, and is clearly reviewed in [START_REF] Cherubini | Application of Probability Distributions to the Modeling of Biogenic CO 2 Fluxes in Life Cycle Assessment[END_REF]. Authors generally agree on the need to account for the complete carbon cycle including biogenic carbon [START_REF] Rabl | How to Account for CO 2 Emissions from Biomass in an LCA[END_REF]Müller-Wenk and Brandão, 2010, etc.), but methodologies vary. Straight-line amortisation of pulse emissions are widely found both in scientific papers [START_REF] Searchinger | Use of U.S. Croplands for Biofuels Increases Greenhouse Gases Through Emissions from Land-Use Change[END_REF]Gnansonou et al., 2009) and methodological guidelines (IPCC, 2006;EU, 2009;PAS2050, 2011, while more dynamic approaches such as time-distributed or delayed emissions have been also developed notably within the LCA framework (Kendall et al., 2009;[START_REF] Müller-Wenk | Climatic Impact of Land Use in LCA-carbon Transfers between Vegetation/soil and Air[END_REF][START_REF] Cherubini | Influence of Allocation Methods on the Environmental Performance of Biorefinery products-A Case Study[END_REF]. As highlighted by some authors, time horizons or time preference in GHG assessments can have a great influence on the final balance [START_REF] Searchinger | Use of U.S. Croplands for Biofuels Increases Greenhouse Gases Through Emissions from Land-Use Change[END_REF][START_REF] Cherubini | Influence of Allocation Methods on the Environmental Performance of Biorefinery products-A Case Study[END_REF]; but they are more related to policy choices rather than a scientific question [START_REF] Fearnside | Time Preference in Global Warming Calculations: A Proposal for a Unified Index[END_REF]. Time horizons already are intrinsically embedded in global warming potentials provided by IPCC (2006) and correspond to a 0.9% annual discount rate, i.e. no preference for time [START_REF] Fearnside | Time Preference in Global Warming Calculations: A Proposal for a Unified Index[END_REF][START_REF] Cherubini | Application of Probability Distributions to the Modeling of Biogenic CO 2 Fluxes in Life Cycle Assessment[END_REF]. We argue that a dynamic accounting for the carbon cycle at the global level would be more scientifically robust but is, for the moment, not in agreement with IPCC methodologies (both the guidelines and global warming potentials). It does neither address the issue of irreversible losses. We also argue that it is a policy issue to allocate the burden of land clearing to the immediate following land use, while driven development factors for a specific land area might be assessed at a more national or even global level. A dynamic accounting for the carbon cost would make PalmGHG rather complex for the users, while producers already had difficulty to adopt the tool due to the policy argument on the above-mentioned development issue related to land clearing. PalmGHG is a management tool that allows for comparing production systems with different management practices. It correctly points to main emission sources and gives incentives to limit land clearing or peat cultivation for instance in order to reduce GHG emissions. Given the mentioned shortcomings, we want to stress that PalmGHG does not provide absolute GHG balances but estimates that are suitable to adapt the practices towards sustainable systems. The RSPO criterion on GHG monitoring comes together with further requirements such as conservation of High Conservation Value areas, which are a more efficient tool to limit irreversible losses such as biodiversity ones.

The idea behind "conservation area" included in the GHG balance is to give further incentive to preserve land area also from a GHG point of view. However, the methodology to measure the carbon savings and include the estimates in the GHG balance has not yet been established. It will face the same issue for biogenic carbon and should be peer-reviewed before being implemented in PalmGHG.

Integrating the spatial and temporal dimensions of the perennial palm crop cycle within a snapshot assessment is not immediate. In PalmGHG, this difficulty is partly overcome by considering data from all oil palm plantings within the whole supply area. Despite the large areas that are often involved, the ages of the oil palms may not be evenly distributed and this may bias the assessment.

Several studies have focused on the treatment and uses of co-products [START_REF] Yacob | Baseline study of methane emission from open digesting tanks of palm oil mill effluent treatment[END_REF][START_REF] Chavalparit | Options for environmental sustainability of the crude palm oil industry in Thailand through enhancement of industrial ecosystems[END_REF][START_REF] Vijaya | Life cycle inventory of the production of crude palm oil -A gate to gate case study of 12 palm oil mills[END_REF][START_REF] Stichnothe | Life cycle assessment of two palm oil production systems[END_REF]. However, the implementation of these technologies involves high capital costs which may restrict their uptake by producers. In addition, there may be limited options for the sale of excess energy especially when grid connection is not possible. However, such technologies can be implemented through clean development mechanisms provided care is taken to avoid double-counting of GHG savings, such as the credits for coal substitution by sale of shell. Moreover, more research effort is needed to better assess fertilising efficiency of land clearance residues and environmental emissions of down-stream processes related to co-product treatment and transport.

The GHG balance is only one of several potential impacts of palm oil production on the environment. However, it has yet to be fully accorded due importance as part of RSPO P&C. In this respect, PalmGHG is a useful tool that can help identify potential GHG savings at the plantation and mill. Together with the other RSPO P&Cs that represent a broader view of sustainability criteria, it can help improve the image of oil palm production by openly accounting for its C footprint and demonstrating improvements leading to greater sustainability.

Compared to other vegetable oils, palm oil and its derivatives usually perform better in terms of GHG balance due to their high yields per hectare. Thus for rapeseed oil [START_REF] Schmidt | Life cycle inventory of rapeseed oil and palm oil[END_REF] Ester respectively. However, comparisons using a single criterion may induce problem shifting in environmental impacts. More comprehensive LCAs must also be undertaken to quantify other environmental impacts and allow for a sounder and more harmonised comparison between the different agricultural products including land transformation and land occupation compared to reference land uses [START_REF] Milà I Canals | Land use impact assessment of margarine[END_REF]. In particular, consideration of impacts on soil fertility and biodiversity is paramount. In this case, other stages of palm oil production might also play an important role such as the use of pesticides causing ecosystem toxicity or boiler emissions affecting human health [START_REF] Schmidt | Life cycle inventory of rapeseed oil and palm oil[END_REF][START_REF] Choo | Determination of GHG contributions by subsystems in the oil palm supply chain using the LCA approach[END_REF][START_REF] Bessou | LCA of Palm Oil in Sumatra, Comparison of Cropping Systems[END_REF].

Conclusions

PalmGHG is a comprehensive GHG calculator representative of the state of the art in terms of available data and international methodologies for GHG accounting. Emphasis has been placed on information directly relevant to palm oil production that should be easily available at the field and mill levels. However, provision is also made for data which might not be available by their substitution with representative default data. Flexibility is an important feature of PalmGHG, with options that allow for alternative calculations and methodology.

During pilot testing it was shown that PalmGHG can identify GHG emission 'hot spots', and so help to define GHG reduction strategies. Feedback from the pilot companies highlighted problems in collecting data, especially those for three consequent years. It should, however, be noted that difficulties related to data recording will progressively diminish once the monitoring of GHG emissions becomes routine. PalmGHG is being reprogrammed to make it more user-friendly in terms of data entry and output display. The new software will allow for the users to quickly generate results and yearly update the database, while being able to readily change default parameters and easily undertake tests of alternative scenarios. On the other hand, difficulties encountered when collecting data for outgrowers are not so easily resolved and indicate a need for a specific strategy to help them record and collect data on a routine basis.

The results of the pilot and scenario tests provide an important information base on which to make recommendations to the RSPO EB and to communicate to a larger audience the work of the RSPO GHG WG and the use of PalmGHG. Further recommendations of the GHG WG to the EB refer to conditions that should be met by new plantations in order to ensure low GHG emissions during their operation.

  PalmGHG approach and boundariesThe PalmGHG calculator is based on a Life Cycle Assessment (LCA) approach; it is hence in line with LCA international standards, i.e. ISO 14040 and 14044. It is more specifically based on an attributional LCA approach(Zamagni et al., 2012) (i.e. the impacts are those linked to the production unit without considering marginal impacts on other productions or any feedback mechanisms), it thus notably excludes indirect land use changes. Indirect land use changes are still not included in most international standards (IPCC, PAS2050, ISCC, GHG Protocol Product Standard, etc.). While LCA norms provide a harmonised global framework to assess environmental impacts over a supply chain, IPCC guidelines provide specific algorithms and emission coefficients to specifically calculate the GHG balance. PalmGHG, like most international GHG accounting standards such as ISCC or PAS2050, is hence based on IPCC(2006). Net GHG emissions due to direct land use change are calculated according toIPCC (2006) and hence account for provisional stocks in the vegetation biomass. Further details on specificities related to perennial cropping systems were provided in PAS2050 by BSI, 2011 (see section 2.2.1).

FiguresFigure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figures

  The results show that high emissions result from clearing logged forest, and conversely that very low (negative) emissions result from clearing low biomass land such as grassland on mineral soils. Fertiliser emissions, including both upstream and field emissions are a not negligible contributor, especially in scenario 3 where net sequestration (sequestration minus land clearing emissions) is high, and in scenarios 4 and 5, where net sequestration is almost null and methane is captured.

	ABOUT HERE
	3.2. Scenario testing
	PalmGHG readily allows manipulation of input data to test management interventions. Results of scenario
	testing are given for a set of dummy data for a base case with the following characteristics: mixed previous
	land uses, 3% of the area containing peat, no POME treatment, OER 20%, own crops mean yield
	20 t FFB/ha, outgrowers' mean yield 14 t FFB/ha (Scenario 1 in Figure 4). The base case scenario (1) repre-
	sents a virtual mill with widespread characteristics for common mills in Indonesia and Malaysia, the first
	world producing countries. It was established by experts who have decades of experience in oil palm agron-
	omy and data collection at field and mill levels. This base case scenario and the tested scenarios are very
	likely to represent existing cases and to address most pregnant interrogations for future plantations (e.g. what
	it the influence of previous land uses on the GHG balance?: scenarios 2-3), for operating mills (e.g. what is

the GHG saving potential through capture and conversion of methane?: scenario 4), and for replanting since main palm plantations in these countries are in their third to fourth generation (e.g. what are the best options to replant with a GHG balance close to zero?: scenario 5).

  reported 5-17 t CO 2 e/t Rapeseed oil and Thamsiriroj and Murphy (2009) and Achten et al. (2010a, b) calculated 39-88 g CO 2 e/MJ Palm Methyl Ester compared to 62 and 124-159 g CO 2 e/MJ for rapeseed and Jatropha Methyl

Table 1 .

 1 Modifications to PalmGHG are still being made, notably to amend default values. Since PalmGHG is meant and robust empirical models and new scientific knowledge. RSPO Emission Reduction Working Group is presently working on establishing specific guidance to check for consistency of user-defined parameters and overall calculations, as well as procedures to program and validate regular updates of PalmGHG. Zamagni, A., Guinée, J., Heijungs, R., Masoni, P., Raggi, A., (2012) Lights and shadows in consequential LCA. Int. J. Life Cycle Assess. 17(7), 904-918. Carbon stocks (below-and above-ground biomass) for the land uses available in PalmGHG Calculated using the OPPRODSIM and OPCABSIM models(Henson, 2005b[START_REF] Henson | Modelling carbon sequestration and greenhouse gas emissions associated with oil palm cultivation and land-use change in Malaysia. A re-evaluation and a computer model[END_REF], with actual amounts dependant on crop age and growth conditions (vigorous or average). The value given here is for a 30 year old, vigorous stand.[START_REF] Harris | Land Use Change and Emission Factors: Updates Since the RFS Proposed Rule[END_REF] pers. com. N. Harris Data 2000 to 2007 Food crops 9 Average of annual and perennial crops in Papua New Guinea where mixed cultivation s practised (WINROCK model, Harris et al., 2009 and Harris pers.

	Tables
	Land uses	Carbon stocks	Reference and origin of data
		in tonne C/ha
	Primary forest	225	Mean of 62 values from the literature adjusted to include roots (Henson, 2012,
			unpublished database)
	Logged forest	87	Henson, 2005a; Henson, 2009
	Coconut	75	European Commission, 2009
	Rubber	62	Yew, 2000; Yew and Nasaruddin, 2002; Henson, 2009, Data are for a mature
			(29 year old) stand and were adjusted to allow for C sequestered in harvested
			rubber wood.
	Cocoa under shade	70	Lasco et al., 2001. Data were adjusted to include roots
	Oil palm	≥ 50
	Secondary regrowth	48	Average of logged forest and food crops
	Shrub	26	WINROCK model (com.)
	Grassland	5	Henson, 2009

to support enlightened management, emphasis was put on including practice-related parameters and enabling field-measurement input data, wherever possible. However, numerous parameters are still not enough sensitive to site-specific conditions (e.g. N 2 O IPCC coefficients, POME or Peat-related emissions, biomass stands etc.). PalmGHG should be updated regularly to introduce updated parameters based on future reported effi-cient

Table 2 .

 2 Pilot mills, their main characteristics and GHG balances assessed with PalmGHG

Mills 7 Mean yield t FFB/ha Outgrowers included Peat soil proportions (own-growers only) Previous land uses t CO 2 e/t CPO

	A	23.6	no	0%	Shrub	-0.02
	B	22.8	no	0%	Shrub	0.01
	C	19.2	no	0%	Logged forest, grassland	0.05
	D	17.2	yes	0%	Mixture of previous land uses	0.45
	E	24.3	no	0%	Cocoa, oil palm	0.83
	F	17.2	yes	0.25%	Logged forest	1.33
	G	15.3	yes	0%	Mixture of previous land uses	1.38
	H	19.3	no	0%	Logged forest, oil palm	1.97
	I	24.2	yes	31%	Grassland, shrub	2.39
	J	18.5	yes	79%	Grassland, shrub	8.32

GWAPP : Global Warming Assessment of Palm Oil Production

PalmGHG may be downloaded under request on the webpage: http://www.rspo.org/en/rspo_palmghg_calculator

Acknowledgements

The authors are very grateful to all who participated in the second RSPO GHG WG for their contributions.

We thank the companies who took part in the pilot phase for supplying data and feedback on the tool. Lastly, we are most grateful to members of the RSPO Executive Board for their continuing support; particular thanks being due to Dr Simon Lord for proposing the adoption of the GWAPP model and to Dr Timothy Killeen for guiding the development of the GHG tool. We also wish to thank the reviewers, whose comments have led to major improvements in the paper.


Daniel Delahaye

Introduction

When an aircraft ies from a city A to a city B, it has to be managed by air traf c controllers in order to avoid collisions with others aircraft. Everyday, about 8000 aircraft y in the French airspace, inducing a huge amount of control workload. Such a workload, is then spread by the mean of the airspace sectoring (the airspace is divided into geometrical sectors, each of them being assigned to a controller team). When ac o n ict between two (or more) aircraft is detected, the controller changes their routes (heading, speed, altitude) in order to keep a minimum distance between them during the crossing. All ying aircraft are then monitored during their navigation and so from departure till destination. The controller has a 2D visualization of his traf c by the mean of the radar screen. At the beginning of radar air traf c control, civil aviation authority used primary radar. In such a system, the radar send a high energy electromagnetic pulse (several mega watts (peak power) during one micro second) in the airspace waiting the echo back. This pulse propagates at the light speed (c = 300 km per second), touch the aircraft and comes back to the radar with a very low power (some nano watts). Based on the time difference between pulses, it is easy to compute the distance between the radar and the aircraft (∆ distance = ∆ time .c 2 where c is the light speed). This system gives a straight measure of the distance between the aircraft and the radar which is a mix between the 2D geographical distance and the altitude of the aircraft. When controllers used only primary radar, they had to asked the pilot to give his altitude and his ID in order to identify more precisely the spot they saw on their screens. The secondary radar really improved this system by the mean of active target concept. In this system, the aircraft is not passive anymore (from the radar point of view) but is equipped by an emitter which answers to the radar when the radar beam reaches the aircraft. Every time the aircraft is addressed by the radar, the aircraft emitter send back two codes called ModeA and ModeC. Those numbers are coded into octal and give respectively the ID and TAS AND WIND ESTIMATION FROM RADAR DATA 2.B.5-1 the altitude of the aircraft. When the radar receives back those two pieces of information, it is easy to compute the 2D geographical distance and the altitude of the aircraft. Another radar system, called Mode S, proposes to establish a data link between the radar and the aircraft during the time the aircraft is in the radar beam. This short connection enables to down link some ight data to the ground such like roll angle, true air speed, heading, altitude rate, etc....and enable to improve the radar tracking [START_REF] Lefas | Improved tracking with mode-s data-linked velocity measurements[END_REF]. This communication between aircraft and the ground can be extended by the mean of satellites which are even accessible over the oceans where there is no radar.

All those radars deliver raw measures with additional noises coming from electromagnetic perturbations, electronic systems, etc.... Those measures are then ltered by the mean of a tracker lter in order to produce smooth trajectories.

When a controller observes its traf c on the radar screen, he tries to identify convergent aircraft which may be in con ict in a near future, in order to apply maneuvers that will separate them. The problem is to estimate where the aircraft will be located in this near future (5-10 minutes); this process is call trajectory prediction. This prediction may be also very useful in order to estimate the workload level in control sector to prevent over capacity event. As a mater of fact, it is very useful to estimate when an aircraft will enter a sector in order to compute the associated sector workload and to apply regulation if necessary. When a sector is expected to be overloaded, the aircraft involved in such a process will be speed up or slow down by the controller in order to adapt the demand to the actual capacity as much as possible. The trajectory prediction depends mainly on the residual noise after ltering (see [START_REF]Con ict detection and resolution in the presence of prediction error[END_REF]): the weight of the aircraft, the temperature and the wind. The residual noise is integrated with time with a growing covariance matrix indicating that the estimated position is less and less accurate. The weight of the aircraft is relevant in the ight dynamic model but is still a raw data. The engines of aircraft are sensitive to the air temperature and such a data is very useful to model the trust of the aircraft b u ti ti sa l s ov e r yd i f cult to measure on real time. Finally, the wind in uences strongly the cinematic of the aircraft and limits also the trajectory prediction. Based on the available accuracy, the actual limit of the trajectory prediction is about 20 minutes. It means that after 20 minutes the uncertainty is so big that the estimated position is no more useful for any ATM purposes.

Several efforts have been tried to improved the trajectory prediction by estimating or suppressing the wind [START_REF] Mondoloni | Improving trajectory forecasting through adaptive ltering technique[END_REF], [START_REF]Wind prediction accuracy for air traf c management decision support tools[END_REF], [START_REF] Rekkas | Three dimensional tracking using on-board measurements[END_REF], [START_REF] Delahaye | Wind eld update using radar track data[END_REF], [START_REF] Bradford | Using aircraft radar tracks to estimate winds aloft[END_REF]. The presents paper, proposes a new method for estimating the wind around aircraft by the mean of observations of the radar tracks and some down linked data. One goal of this work is to show how it is possible to extract wind information from the radar observations.

The paper is organized as follow : the rst part presents the relation between air speed vector, ground speed vector and the wind. The second part gives the observability condition of the wind based on the available measures. The third part present the Kalman lter and its non linear extensions. The fourth part presents the different models which have been used for the wind extraction. The fth part presents some results and compares the performance of our models. Finally, the sixth part describes how a wind eld map can be built using some spread measures of the wind.

General Relations between Speed Vectors

The following notations will be used in the paper. Vectors and matrices are shown with underlined symbols. Speeds

V = v x v y V = v 2 x + v 2 y θ g = arctan v x v y
Where θ g is the The route angle (with reference to north). The same notation will be used for the air speed (T ,t x ,t y , θ a ) and the wind (W , w x , w y , θ w ). We have :V = T + W . These relations are shown on gure 1. dt . One can show that the air turning rate (ω a ) and the ground turning rate (ω g )a r e related by the following expression :

Air Turning Rate

ω a = dθ a (t) dt Ground Turning Rate ω g = dθ g (t) dt = d[arctan u(t)] dt = d[arctan u] du . du dt with u = v x v y ω g = γ x .v y -v x .γ y V 2 (1)
ω g = T 2 +T * W * cos(θ a (t)-θ g (t)) T 2 +W 2 +2.T.W cos(θ a (t)-θ g (t)) ω a (2)
Hypothesis We consider en-route traf c( t r a f c away from airports), with aircraft ying in cruise phase. Aircraft are supposed to y at constant air speed (T = C te ) and turn with constant air turning rate (ω a = C te ). The average wind is supposed to be constant in the neighborhood of the aircraft.

Observability Conditions

Mode S Radar

Like it has been mentioned in the introduction, such radar are able to access on board parameters and especially airspeed vector. Having such measures the system is fully determined at any time as it is shown in the following equations. (3)

In this system there are two equations and two unknowns (w x , w y ); the other parameters are given by the radar (v x , v y , T, θ a ). A Kalman lter will be used to extract the wind by removing the noise.

Standard Radar

If measures come from classical radars, only position measures are available (x, y) and two situations have to be taken into account.

First situation: One aircraft

In the case where radar measures come from one aircraft, wind may be observed only after two asymmetric turns (meaning three straight lines separated by two turns). As a mater of fact, the rst segment brings the following system of equations ; v 1x = Tsin(θ a1 )+w x v 1y = Tcos(θ a1 )+w y [START_REF] Bradford | Using aircraft radar tracks to estimate winds aloft[END_REF] In this system (4), there are 4 unknowns (T, θ a1 , w x , w y ) and two equations, so two equations are missing.

After the rst turn two new equations are added to this system with one extra unknown (θ a2 ):

v 1x = Tsin(θ a1 )+w x v 1y = Tcos(θ a1 )+w y v 2x = Tsin(θ a2 )+w x v 2y = Tcos(θ a2 )+w y ( 5
)
In this new system (5), there are 5 unknowns (T, θ a1 , θ a2 , w x , w y ) and only 4 equations. It is only after the second turn that the system is fully determined (6 unknowns (T, θ a1 , θ a2 , θ a3 , w x , w y ) and 6 equations) :

⎧ ⎨ ⎩ v 1x = Tsin(θ a1 )+w x v 1y = Tcos(θ a1 )+w y v 2x = Tsin(θ a2 )+w x v 2y = Tcos(θ a2 )+w y v 3x = Tsin(θ a3
)+w x v 3y = Tcos(θ a3 )+w y (6) This system has a closed form solution for which the wind is given by :

w x = (v 3y -v 2y )V 2 1 +(v 1y -v 3y )V 2 2 +(v 2y -v 1y )V 2 3 2{v 1y (v 2x -v 3x )+v 2y (v 3x -v 1x )+v 3y (v 1x -v 2x )} (7)
w y = (v 2x -v 3x )V 2 1 +(v 3x -v 1x )V 2 2 +(v 1x -v 2x )V 2 3 2{v 1y (v 2x -v 3x )+v 2y (v 3x -v 1x )+v 3y (v 1x -v 2x )} (8)
With

V 1 = v 2 1x + v 2 1y V 2 = v 2 2x + v 2 2y V 3 = v 2 3x + v 2 3y
For both expressions (w x , w y ) the denominators must not be equal to zero meaning that turns have to be asymmetric 2.B.5-3

(V 1 = V 2 = V 3 ).

Figure 1. Speed Vector Relations

It is also possible to extract a close form for the True Air Speed (T ) of the aircraft using only the ground speeds :

T = ∆ v 12 . ∆ v 13 . ∆ v 23 2.|v 1y .(v 2x -v 3x )+v 2y .v 3x -v 2x .v 3y +v 1x .(v 3y -v 2y )| (9)
Where ∆ v ij = v iv j T h es a m ek i n do fc l o s e forms can be extracted for the three heading angles (θ a1 , θ a2 , θ a3 ):θ ai = arctan t ix t iy where

t ix t iy = 2.v ix .A + B 2.v iy .A -C
where

A = ∆ v 23x .v 1y + ∆ v 31x .v 2y + ∆ v 12x .v 3y B = ∆ v 23y .V 2 1 + ∆ v 31y .V 2 2 + ∆ v 12y .V 2 3 C = ∆ v 23x .V 2 1 + ∆ v 31x .V 2 2 + ∆ v 12x .V 2 3
Second situation : Two aircraft When radar measures are available for two aircraft (a and b), only one turn for both trajectories is needed to have enough information for wind estimation. The two rst segments bring the following system of equations with 6 unknowns (T a , T b , θ a a1 , θ a b1 , w x , w y ) and 4 equations : va 1x = T a sin(θ a a1 )+w x va 1y = T a cos(θ a a1 )+w y vb 1x = T b sin(θ a b1 )+w x vb 1y = T b cos(θ a b1 )+w y (10) After the second turn, the new systems is fully determined with 8 unknowns(T a , T b , θ a a1 , θ a b1 , θ a a2 , θ a b2 , w x , w y )a n d8 equations :

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ va 1x = T a sin(θ a a1 )+w x va 1y = T a cos(θ a a1 )+w y vb 1x = T b sin(θ a b1 )+w x vb 1y = T b cos(θ a b1 )+w y va 2x = T a sin(θ a a2 )+w x va 2y = T a cos(θ a a2 )+w y vb 2x = T b sin(θ a b2 )+w x vb 2y = T b cos(θ a b2 )+w y (11)
The associated closed form of the wind is given by :

w x = (vb 1y -vb 2y )(Va 2 1 -Va 2 2 )+(va 2y -va 1y )(Vb 2 1 -Vb 2 2 ) 2{(va 1x -va 2x )(vb 1y -vb 2y )-(va 1y -va 2y )(vb 1x -vb 2x )} (12) w y = (vb 2x -vb 1x )(Va 2 1 -Va 2 2 )+(va 1x -va 2x )(Vb 2 1 -Vb 2 2 ) 2{(va 1x -va 2x )(vb 1y -vb 2y )-(va 1y -va 2y )(vb 1x -vb 2x )} (13)
with

Va 1 = va 2 1x + va 2 1y Va 2 = va 2 2x + va 2 2y Vb 1 = vb 2 1x + vb 2 1y Vb 2 = vb 2 2x + vb 2 2y and (V a 1 = V a 2 ;V b 1 = V b 2 ).
Like for the one aircraft case, it is very easy to extract the other unknowns :

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ta 1x = va 1x -w x ta 1y = va 1y -w y tb 1x = vb 1x -w x tb 1y = vb 1y -w y ta 2x = va 2x -w x ta 2y = va 2y -w y tb 2x = vb 2x -w x tb 2y = vb 2y -w y then ⎧ ⎨ ⎩ T a = ta 2 1x + ta 2 1y = ta 2 2x + ta 2 2y T b = tb 2 1x + tb 2 1y = tb 2 2x + tb 2 2y
Ac l o s ef o r mo fT a and T b can be extracted : Based on those observability conditions, several wind estimation models have been developed for which some of them use a Kalman lter. This lter is now described in the following section.

T a = ta 2 1x + ta 2 1y = ta 2 2x + ta 2 2y ta 2 1x = ∆vb 12y Va 2 1 -2.∆va 12x .va 1x -Va 2 2 +

Kalman Filtering

Linear Form

The Kalman lter is named after Rudolph E.Kalman, who in 1960 published his famous paper describing a recursive solution to the discrete-data linear ltering problem [START_REF] Kalman | A new approach to linear ltering and prediction problems[END_REF], [START_REF] Grewal | Kalman Filtering Theory and Practice[END_REF], [START_REF] Welch | An introduction to the kalman lter[END_REF]. The Kalman lter is essentially a set of mathematical equations that implements a predictor-corrector type estimator that is 2.B.5-4 optimal in the sense that it minimizes the estimated error covariance (when some presumed conditions are met). It has been successfully applied in many real problems [START_REF] Maybeck | Stochastic Models, Estimation, and Control[END_REF], [START_REF] Chui | Kalman Filtering with Real-Time Applications[END_REF], [START_REF] Battin | Astraunotical Guidance[END_REF], [START_REF] Asher | Filtering for precision pointing and tracking with application for aircraft to satellite tracking[END_REF], [START_REF] Anderson | Optimal Filtering[END_REF]. The following description of the Kalman lter is mainly due to G.Welch and G.Bishop [START_REF] Welch | An introduction to the kalman lter[END_REF].

The Kalman lter addresses the general problem of trying to estimate the state X (k) ∈ R n of a discrete controlled process that is governed by the linear stochastic difference equation call evolution equation :

X (k + 1)=F(k).X(k)+G(k). Û(k)+v(k) ( 14
)
with a measurement Z(k) ∈ R m :

Z(k)=H(k).X(k)+w(k) ( 15
)
The n.n matrix F(k) in the evolution equation [START_REF] Maybeck | Stochastic Models, Estimation, and Control[END_REF] relates the state at the present step k to the state at the next step k +1, in the absence of either a driving function or process noise. The n.p matrix G(k) relates the optional control input Û(k) ∈ R p to the state X .T h e n.m matrix H(k) in the measurement equation [START_REF] Mondoloni | Improving trajectory forecasting through adaptive ltering technique[END_REF] relates the state to the measurement Z(k). The random variables v(k) and w(k) represent the process and measurement noise (respectively). They are assumed to be independent (of each other), centered, white, and with normal probability distributions. Furthermore the optional control is disturbed by a random variable n(k) ( Û(k)=U(k)+n(k)) which will be supposed to have a normal probability distribution also.

⎧ ⎨ ⎩ p(v(k)) ∽ N (0, Q(k)) p(w(k)) ∽ N (0, R(k)) p(n(k)) ∽ N (0, N(k)) (16)
So we have :

E [v(k)] = 0 E v(k)v( j) T = Q(k).δ kj E [w(k)] = 0 E w(k)w( j) T = R(k).δ kj E [n(k)] = 0 E n(k)n( j) T = N(k).δ kj ( 17
)
where T is the transpose operator and δ kj is the Kronecker symbol which is 1 if k = j a n d0i fk = j.

We de ne X(k + 1/k) ∈ R n to be the a priori state estimate at step k + 1 given knowledge of the process prior to step k + 1, and X(k + 1/k + 1) ∈ R n to be the a posteriori state estimate at step k + 1g i v e n measurement Z(k + 1). We can then de ne a priori and a posteriori estimate error as :

e(k + 1/k)=X (k + 1) -X(k + 1/k) e(k + 1/k + 1)=X (k + 1) -X(k + 1/k + 1) (18)
The a priori estimate error covariance is then :

P(k + 1/k)=E e(k + 1/k)e(k + 1/k) T (19)
and the a posteriori estimate error covariance is :

P(k + 1/k + 1)= E e(k + 1/k + 1)e(k + 1/k + 1) T (20)
In deriving the equation for the Kalman lter, we begin with the goal of nding an equation that computes an a posteriori state estimate X(k +1/k +1) as a linear combination of an a priori estimate X(k + 1/k) and a weighted difference between an actual measurement Z(k +1) and a measurement prediction H. X (k +1/k) :

X(k + 1/k + 1)= X (k + 1/k)+ K(k + 1). Z(k + 1) -H. X (k + 1/k) (21)
The difference Z(k + 1) -H. X (k + 1/k) in equation 21 is called the measurement innovation, or the residual. The residual re ects the discrepancy between the predicted measurement H. X (k +1/k) and the actual measurement Z(k + 1). A residual of zero means that the two are in complete agreement.

The n.m matrix K(k + 1) in equation 21 is chosen to be the gain or blending factor that minimizes the a posteriori error covariance 20. This minimization can be accomplished by rst substituting equation 21 into the above de nition for e(k + 1/k + 1), substituting that into equation 20, performing the indicated expectations, taking the derivative of the trace of the result with respect to K(k + 1), setting that result equal to zero, and then solving for K(k + 1). One form of the resulting K(k +1) that minimizes equation 20 is given by :

K(k + 1)=P(k + 1/k).H(k + 1) T . H(k + 1)P(k + 1/k)H(k + 1) T + R(k + 1) -1 (22)
Looking at equation 22 we see that as the measurement error covariance R(k + 1) approaches zero, the gain K(k + 1) weights the residual more heavily. Speci cally : lim

R(k+1)→0 K(k + 1)=H -1 (23)
On the other hand, as the a priori estimate error covariance P(k +1/k) approaches zero, the gain K(k +1) weights the residual less heavily. Speci cally : lim

P(k+1/k)→0 K(k + 1)=0 (24)
2.B.5-5

As the measurement error covariance R(k + 1),t h e actual measurement Z(k + 1) is "trusted" more and more, while the predicted measurement H. X (k + 1/k) is trusted less and less. On the other hand, as the a priori estimate error covariance P(k +1/k) approaches zero the actual measurement Z(k + 1) is trusted less and less, while the predicted measurement H. X(k + 1/k) is trusted more and more. The Kalman lter estimates a process by using a form of feedback control : the lter estimates the process state at some time and then obtains feedback in the form of (noisy) measurements. As such, the equation for the Kalman lter fall into two groups : time update equations and measurement update equations. The time update equations are responsible for projecting forward (in time) the current state and error covariance estimates to obtain the a priori estimates for the next time step. The measurement update equations are responsible for the feedback -i.e. for incorporating a new measurement into the a priori estimate to obtain an improved estimate.

The time update equations are the following :

⎧ ⎨ ⎩ X(k + 1/k)=F(k).X(k)+G(k).U(k) P(k + 1/k)=F(k).P(k/k).F(k) T + G(k).N(k).G(k) T + Q(k) (25)
The rst equation is the state prediction and the second one gives the associated covariance. The measurement update equations are the following :

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ K(k + 1)=P(k + 1/k).H(k + 1) T . H(k + 1)P(k + 1/k)H(k + 1) T + R(k + 1) -1 X (k + 1/k + 1)= X(k + 1/k)+K(k + 1). Z(k + 1) -H(k + 1). X(k + 1/k) P(k + 1/k + 1)=[I -K(k + 1).H(k + 1)]. P(k + 1/k) (26)
The rst task during the measurement update is to compute the Kalman gain, K(k + 1). The next step is to actually measure the process to obtain Z(k + 1), and then to generate an a posteriori state estimated by incorporating the measurement. The nal step is to obtain an a posteriori error covariance estimate.

After each time and measurement update pair, the process is repeated with the previous a posteriori estimates used to project or predict the new a priori estimate. In order to initiate the process, the rst initial state and the associated covariance has to be set up. Usually, the rst state is build with an initial measure and the associated covariance is set up with the measure covariance or with a big values, because the lter has no con dence in its initial state (it comes from a measure). This Kalman lter has been used in our wind estimation problem for which the associated models are now presented.

Extended Form

The standard Kalman lter addresses the general problem of trying to estimate the state X (k) ∈ R n of a discrete-time controlled process that is governed by a linear stochastic difference equation. But what happens if the process to be estimated and (or) the measurement relationship to the process is non-linear ? A Kalman lter that linearized about the current mean and covariance is referred to as an extended Kalman lter or EKF. We can linearize the estimation around the current estimate using the partial derivative of the process and measurement functions to compute estimates even in the face of non-linear relationships. Let us assume that our process again has a state vector X(k) ∈ R n , but that process is now governed by the non-linear stochastic difference equation :

X(k + 1)=F [k, X(k),U(k)] + v(k) (27)
with a measurement Z(k) ∈ R m that is :

Z(k + 1)=H [k, X(k)] + w(k) ( 28
)
where the random variable v(k) and w(k) again represent the process and measurement noise. The new time update equations are now the followings :

⎧ ⎨ ⎩ X(k + 1/k)=F [k, X(k),U(k)] P(k + 1/k)=F X (k).P(k/k).F X (k) T + F U (k).N(k).F X (k) T + Q(k) (29)
where F X (k) is the Jacobian matrix of partial deriva- tives of F with respect to the vector X and is given by the following :

F X (k) = ▽ X F [k, X,U] T T X= X(k/k),U= Û(k) (30)
where

▽ S = ∂ ∂ S 1 , ∂ ∂ S 2 ,..., ∂ ∂ S n s (31)
is the differential operator.

2.B.5-6

The linearized measurement update is now given by :

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ K(k + 1)=P(k + 1/k).H X (k + 1) T . H X (k + 1).P(k + 1/k)H X (k + 1) T + R(k + 1) -1 X(k + 1/k + 1)= X (k + 1/k)+K(k + 1). Z(k + 1) -H k, X(k + 1/k) P(k + 1/k + 1)= I -K(k + 1).H X (k + 1) . P(k + 1/k)
(32) Unlike its linear counterpart, the extended Kalman lter is not an optimal estimator. In addition, if the initial estimate of the state is wrong, or if the process is modeled incorrectly, the lter may quickly diverge, owing to its linearization. Another problem with the extended Kalman lter is that the estimated covariance matrix tends to underestimate the true covariance matrix and therefore risks becoming inconsistent in the statistical sense without the addition of "stabilizing noise". When the state transition and observation models are highly non-linear, the extended Kalman lter can give particularly poor performance [START_REF] Julier | A new extension of the kalman lter to nonlinear systems[END_REF].This is because only the mean is propagated through the non-linearity.

The Unscented Kalman Filter(UKF) is another non linear extension of the linear Kalman lter for which both mean vector and covariance matrix are propagated through the non-linearity. A detailed description of the UKF is given in appendix[A].

Having described our problem and some estimation techniques, the next section introduces the models which haves been used for our experiments.

Models

Model with Air Speed and Air Turning Rate (Model 1)

The state vector is given by :

X(k)=[x(k) y(k)t x (k)t y (k) w x (k) w y (k)] T (33)
where x(k), y(k) is the position, t x (k),t y (k) the True Air Speed (TAS) and w x (k), w y (k) the wind. The measure vector consists in the radar position and the true air speed :

Z(k)=[x m (k) y m (k)t x m (k)t y m (k)] T (34)
Having access to the air turning rate (ω a ), it can be included in the prediction matrix. The structure of the system is the following :

X(k + 1)=F(k).X(k)+v(k) (35)
where

F(k)= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 10 C 1 (ω a ) C 2 (ω a ) ∆ t 0 01-C 2 (ω a ) C 1 (ω a ) 0 ∆ t 00 C 3 (ω a ) C 4 (ω a ) 00 00-C 4 (ω a ) C 3 (ω a ) 00 00 0 0 1 0 00 0 0 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (36) Z(k)= ⎡ ⎢ ⎢ ⎣ 100000 010000 001000 000100 ⎤ ⎥ ⎥ ⎦ .X (k)+w(k) (37)
where

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ C 1 (ω a (k)) = sin(ω a (k)∆ t ) ω a (k) C 2 (ω a (k)) = 1-cos(ω a (k)∆ t ) ω a (k) C 3 (ω a (k)) = cos(ω a (k)∆ t ) C 4 (ω a (k)) = sin(ω a (k)∆ t ) (38)
The model being exact, the model noise covariance matrix R = 0. The measure noise covariance matrix R is given by :

R = ⎡ ⎢ ⎢ ⎣ σ 2 p 000 0 σ 2 p 00 00σ 2 T 0 000σ 2 T ⎤ ⎥ ⎥ ⎦
with σ p = 100 meters and σ T = 0.2k t s .

Model with Air Speed only (model 2)

The state vector is given by :

X(k)=[x(k) y(k)t x (k)t y (k) w x (k) w y (k)] T (39)
The measure vector consists in the radar position and the true air speed :

Z(k)=[x m (k) y m (k)t x m (k)t y m (k)] T (40)
T h ep r e d i c t i o nm a t r i xi sn o wg i v e nb y:

F(k)= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 10∆ t 0 ∆ t 0 01 0 ∆ t 0 ∆ t 0 01000 0 00100 0 00010 0 00001 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (41)
2.B.5-7

The measure matrix H and the associated covariance matrix R are the same as in the rst model. This model is linear but is false for the air speed vector evolution. In order to take into account this model error, the following covariance matrix (Q) is included in the lter : Q is given by :

Q = ⎡ ⎢ ⎢ ⎣ 00 0 0 00 0 0 001 .00 0 001 .0 ⎤ ⎥ ⎥ ⎦

Model with Turning Rate Only (Model 3)

This model is the same as the rst one with the following measure vector ;

Z(k)=[x m (k) y m (k)] T (42)
and the associated measure equation

Z(k)= 100000 010000 .X(k)+w(k) (43)

Models Without Command

Model using Kalman Filter (model 4)

Having no access to the turning rate, this parameter may be extracted by a Kalman lter (ω a included in the state vector). If the turning rate is included in the state vector, the evolution of the system is not linear and an EKF or an UKF has to be used to manage such state vector.

Model using turns (model 5)

As it has been shown in the previous section (observability conditions), wind may be extracted by observing ground radar track during turns.

The key element of this approach is the turning rate detector based on equation [START_REF] Anderson | Optimal Filtering[END_REF]. When ω g is greater than a given threshold, the aircraft is consider to be in turn. Based on this turn detector, straight line segment are easily identi ed for which ground speed vector averages are computed. Speed vector estimates in straight lines are given by the framework of Figure 2. In this gure, the turns are rst detected and a counter is then used to select the right averaging process in order to compute ground speed estimates. When two aircraft are used for wind estimation, the same framework is duplicated and counters select only two averaging blocks in order to produce the following estimates :V a1 ,V a2 ,V b1 ,V b2 .

Results

Simulation Framework

In order to test and compare our models, a radar tracker simulator has been used to produce the reference trajectory which has been disturbed with a Gaussian noise. The disturbed trajectory is then ltered by a Kalman lters which generates the estimated trajectory which is then compared to the reference trajectory. This framework is summarized on Figure 3. For all experiments, a wind of 40 kts has been used with

θ w = 240 o ⇒ W =[-34.64kts, -20kts] T = [-17.82m/s, -10.28m/s] T .

Models with Mode-S radar data (model 1,2 and 3)

The rst trajectory used for our experiments is built with 3 straight lines (20 minutes for each) connected with turns as it can be seen on gure 4.This trajectory has been disturbed by a Gaussian noises for which the means are zero and the standard deviations are the 2.B.5-8 Results for models 1 and 2 Models 1 and 2 have been tried on this trajectory for which the residual errors on the wind estimates are given on gure 5 , 6, 7 and 8. The gures 5 and 6 show the residual error for the wind strength estimate. Figure 5 shows the convergence phase of the lter ( rst 6 minutes) and gure 6 shows the residual error after the convergence phase. The same kind of results are given on gures 7 and 8 for the residual wind angle estimates. The models 1 and 2 have the same evolution during the convergence phase which After this convergence phase, the wind strength estimate error stays below 0.2kts for both models 6. The model 1 has a more accurate behavior at the end of the simulation with an error which stay below 0.1kts. The convergence phase for the wind angle estimate last longer than for the wind strength (about 6 minutes (see gure 7)). After this convergence phase the wind angle residual error stays below 1 degree (in absolute value) for both model. As for the wind strength estimation, the second model (with turning rate) produce a better estimate. The rst model is disturbed by the turns. Both models model are able to produce accurate wind estimate all along the trajectory (after the convergence phase) even during the turns.

Results for models 3 is not able to converge and model 3 (which is models 2 without air speed measure) has to wait the second turn to be able to produce an estimate of the wind (see gure 9 , 10 ,11 and 12). During the convergence phase ( gure 9), model 3 reduces the wind strength error till 40kts which is the limit for such estimation because the Kalman lter is not able to distinguish in the position measures the part coming from the air speed and the part coming from the wind. After the convergence phase (see gure 10) the lter has to wait the second turn in order to be able to produce a reliable estimate of the wind strength (as it has been previously shown (observability conditions), it is only after the second turn that the wind observability conditions are met). As for the wind strength error, the wind angle error reach zero after the second turn.

Models with Classical Radar DataModel 4 and 5)

Results for models 4

When on board measures (TAS, heading, turning rate) are not available, the turning rate has to be included in the state vector inducing a non linear prediction equation of the lter. In order to address such non linear evolution, an EKF and an UKF has been and tested on the rst trajectory. Both lters are not able to estimate the wind vector even after the second turn. This behavior is due to the residual error of the Taylor expansion used in the EKF or in the UKF. Based on the observation (positions) the lter is not able to nd the part due to the wind and the one due to the error coming from the Taylor expansion.

Results for model 5

Based on the turning rate estimate given by equation 1, it is very easy to build a turn detector when |ω g | > threshold. The ground vector speed has been disturbed by a Gaussian noise with zero mean and 0.2kts standard deviation (performance of the actual trackers). This method has been used on the rst trajectory. The ground turning rate estimate is given by gure 13. The two turns are well identi ed with a value of about +/-1.3 degree/sec in turns (it must be noticed that the ground turning rate is different from the air turning rate). The results of this approach are Those results has been computed by using equations [START_REF]Wind prediction accuracy for air traf c management decision support tools[END_REF], [START_REF] Delahaye | Wind eld update using radar track data[END_REF] and the following ground speed estimates coming from the averaging process : Having now some wind estimates on some points in the airspace (where the aircraft are located), the next step consists in the global wind eld interpolation based on the meteorological model called "Shallow-Water". This model is described in [START_REF] Cohn | The behavior of the forcast error covariances for a kalman lter in two dimensions[END_REF] and the associated discret expansion in given in [START_REF] Richtmyer | Difference Methods for Initial-Value Problems[END_REF].

v 1x =

Vector Spline Interpolation

Shallow Water Wind Model

Let V be the vector eld giving the wind velocity at each point. We will assume in the following that area of interest is located high enough above the ground level so that viscous effects are negligible. Assume furthermore that wind eld is two dimensional (that is does not depends on altitude). Then the evolution of V is described by the system of equations : with :

⎧ ⎪ ⎨ ⎪ ⎩
V = U + u v
U the mean west-east component of the wind, f is the Coriolis parameter which can be expressed as :

f = 2Ω sin φ 0 -2ΩR -1 cosφ 0
with Ω = 2π 86400 the rotation rate of the earth, φ 0 the latitude (assumed to be constant in the domain of interest) and R = 6.370 × 10 6 m the earth radius. Φ is the geopotential of the isobaric surface and is written as : Φ = θ(y)+φ with :

θ(y)=θ 0 -2Ω(sinφ 0 )y -ΩR -1 (cosφ 0 )y 2 U
and :

φ 0 = 3 × 10 4 ms -2 U = 25ms -1
The set of equations 44 is known as shallow water model. With little rewriting, it can be shown that this model has the following synthetic form :

( ∂ ∂t + L)X = 0 with : X = ⎡ ⎣ u v φ ⎤ ⎦
and L the differential operator :

L = A ∂ ∂X + B ∂ ∂y +C where A = ⎡ ⎣ U 00 0 U 0 Φ 0 U ⎤ ⎦ B = ⎡ ⎣ 000 001 0 Φ 0 ⎤ ⎦ C = ⎡ ⎣ 0 -f 0 f 00 0 Φ ′ 0 ⎤ ⎦ Interpolating Vector Splines 2.B.5-12
The symbol σ(y, ξ) of this differential operator L easily computes as :

σ(y, ξ)= ⎡ ⎣ -iξ 1 U -f 0 f -iξ 1 U -iξ 2 -iξ 1 Φ(y) -iξ 2 Φ(y)+Φ ′ (y) -iξ 1 U ⎤ ⎦
Since a solution 0 of the Shallow-Water model must ful ll :

( ∂ ∂t + L)V = 0
we will seek for a vector eld V such that :

V = argmin Z T 0 Z D ( ∂ ∂t + L)X (x,t) 2 dxdt
under the constraints that :

X (x i ,t i )=v i , i = 1 ...n
if v i is the wind measurement at time t i and location x i , D is the domain of interest and I is the time interval chosen for computation. Assuming that boundary conditions vanish, we obtain that the functional criterion has expression :

Z I Z D (- ∂ 2 ∂t 2 + L t L)X (x,
σ(y, ξ)= ⎡ ⎣ ξ 2 1 (Φ 2 + U 2 )+ f 2 ξ 1 ξ 2 Φ + iξ 1 (ΦΦ ′ -2 fU) -iξ 1 + ξ 2 U Φ) ξ 1 ξ 2 Φ -iξ 1 (ΦΦ ′ -2 fU f 2 + ξ 2 1 U 2 + ξ 2 2 Φ 2 + Φ ′2 ξ 1 ξ 2 (U + Φ) -iξ 2 Φ ′ iξ 1 + ξ 2 U Φξ 1 ξ 2 (U + Φ)+iξ 2 Φ ′ ξ 2 1 ⎤ ⎦
L t L admits an eigenvector expansion (φ n ) n∈N with associated eigenvalues (λ n ) n∈N . We will week for a kernel K(x, y,t, u) such that : with : n = √ λ n . The constraints that X (x i ,t i )=v i , i = 1 ...n may thus be rewritten as :

Z I Z D K(x i , y,t i , u), (- ∂ 2 ∂t 2 + L t L)X (y, u) dydu = v i
which is much more tractable than the point version. Introducing Lagrange multipliers (ε i ) i=1...n , the problem may be solved and the solution written as :

X (x,t)= ∑ i=1...n ε i K(x i , x,t i ,t) so that optimal vector eld is a linear combination of kernels evaluated at measurement points (x i ,t i ) i=1...n . This is a classical result within the frame of interpolating spline theory, with the restriction that the kernel is not translation invariant. From now, only numeric evaluations of the kernel can be made, based on nding nite elements approximations of the eigenfunctions, then truncating the in nite sum in the kernel expansion at an order where residual is low enough. This yields a piecewise polynomial approximate kernel, which can be computed once and used afterward (note that if the kernel is computable, nding the wind eld expansion is done by solving a linear system).

Conclusion

This paper has given a new approach for extracting the wind information from the radar tracks. Two approaches has been presented. When True Airspeed measures are available, linear models may be used with regular Kalman lter. In this rst approach, wind estimate are available all along the aircraft trajectory. When only position measures are available, an observability analysis has shown that wind may be estimated only after turns (one turn for two aircraft or two turns for one aircraft). For such approach, a closed form of the wind has been developed. Then, those models have been validated with realistic simulations. Having wind measures spread in the airspace, an interpolating method using vector splines has been developed. Based on the differential operator of the Sallow-Water model and the local wind observations, this method build a vector eld using kernel functions. The associated eigenfunctions (φ n )h a v et ob e computed numerically on each point of the grid where the vector eld has to be computed. It must be noticed that the Shallow-Water is only valid above 2000m; this 2.B.5-13 means this method may be used mainly for En-Route traffic.

  v x = Tsin(θ a )+w x v y = Tcos(θ a )+w y

Figure 3 Figure 4 .Figure 5 .

 345 Figure 3. Test Framework

Figure 6 .Figure 7 .

 67 Figure 6. Residual Wind Strength error for

Model 1 and 2 ModelFigure 8 .Figure 9

 289 Figure 8. Residual Wind Angle Error for Model 1 and 2

3 Figure 10 .Figure 11 . 3 Figure 12 .2Figure 13 .

 3101131213 Figure 10. Residual Wind Strength Error for Model 3

Figure 14 .

 14 Figure 14. Second trajectory set. Each straight lines l ast 20 minutes. The turning rate in the rst trajectory is 1 deg and -1 deg for the second one. The true air spee of the rst aircraft is 300kts and 400 for the second one

  x, y,t, u), (-∂ 2 ∂t 2 + L t L)X (y, u) dydu = X (x,t)Formally, we may write :K(x, y,t, u)= ∑ n c n (t, u)φ n (x) ⊗ φ n (y)and nd the right expression for c n (t, u). In the case I = R and wind elds of nite energy over time, some computations show that the sought after kernel has expression :K(x, y,t, u)= ∑ n e -n |t-u|2 n φ n (x) ⊗ φ n (y)

  54.4818 m/s v 1y = 61.9523 m/s v 2x = 84.3536 m/s v 2y = -10.2142 m/s v 3x = -17.6780 m/s v 3y = 91.8504 m/s When two aircraft trajectories are available, only one turn for each is necessary. For such a model, the trajectories given on gure 14 has been used. The turning rate estimates are given on gures 15 and 16. va 1x = 90.5494 m/s va 1y = 98.0082 m/s va 2x = 90.5552 m/s va 2y = -118.4478 m/s vb 1x = -221.7254 m/s vb 1y = -10.2111 m/s vb 2x = -17.6796 m/s vb 2y = -214.2995 m/s

				First aircraft turning rate estimate
		1.2
						Ground turning Rate
		1
		0.8
	Turning Rate (degres/sec)	0.4 0.6
		0.2
		0
		-0.2
		0	500	1000	1500	2000	2500
				Time (sec)
	Figure 15. Turning Rate Estimate of the First Aircraft
			First Trajectory
				Second aircraft turning rate estimate
		0.1
						Ground turning Rate
		0
		-0.1
	Turning Rate (degres/sec)	-0.6 -0.5 -0.4 -0.3 -0.2
		-0.7
		-0.8
		-0.9
		-1
		0	500	1000	1500	2000	2500
				Time (sec)
	Figure 16. Turning Rate Estimate of the Second Aircraft
							The results of this approach are given in the following
							table :
							Wind	Strength	Error
							Est	Angle
							39.64kts w y = -10.2227m/s 239.918deg 0.082deg 0.36kts w x = -17.6485m/s Those results has been computed by the using equa-
							tions 12 13 and the ground speed estimates coming
							from the averaging process :
							2.B.5-11

Appendix A : Unscented Kalman Filter

The unscented Kalman lter (UKF) [START_REF] Julier | A new extension of the kalman lter to nonlinear systems[END_REF] uses a deterministic sampling technique known as the unscented transform to pick a minimal set of sample points (called sigma points) around the mean. These sigma points are then propagated through the nonlinear functions and the covariance of the estimate is then recovered. The result is a lter which more accurately captures the true mean and covariance. In addition, this technique removes the requirement to analytically calculate Jacobians, which for complex functions can be a dif cult task in itself.

Prediction

As with the EKF, the UKF prediction can be used independently from the UKF update.

The estimated state and covariance are augmented with the mean and covariance of the process noise.

A set of 2L + 1 sigma points is derived from the augmented state and covariance where L is the dimension of the augmented state.

where is the ith column of the matrix square root of using the de nition: square root A of matrix B satis es

The matrix square root should be calculated using numerically ef cient and stable methods such as the Cholesky decomposition. The sigma points are propagated through the transition function .

The weighted sigma points are recombined to produce the predicted state and covariance.

where the weights for the state and covariance are given by: Typical values for , and are 10 3 , 2 and 0 respectively.

Update

The predicted state and covariance are augmented as before, except now with the mean and covariance of the measurement noise.

As before, a set of 2L + 1 sigma points is derived from the augmented state and covariance where L is the dimension of the augmented state.

Alternatively if the UKF prediction has been used the sigma points themselves can be augmented along the following lines The sigma points are projected through the observation function

The weighted sigma points are recombined to produce the predicted measurement and predicted measurement covariance.

The state-measurement cross-covariance matrix, is used to compute the UKF Kalman gain.

As with the Kalman lter, the updated state is the predicted state plus the innovation weighted by the Kalman gain, And the updated covariance is the predicted covariance, minus the predicted measurement covariance, weighted by the Kalman gain.