Hakim Bouadi
email: hakimbouadi@yahoo.fr.

Félix Mora-Camino
email: moracamino@hotmail.fr.

Daniel Choukroun
email: d.choukroun@tudelft.nl.

Space-Indexed Control for Aircraft Vertical Guidance with Time Constraint

With the growth of civil aviation traffic capacity, safety and environmental considerations urge today for the development of guidance systems with improved accuracy for spatial and temporal trajectory tracking. This should induce increased practical capacity by allowing timely operations at minimum separation standards while at takeoff and landing, trajectory dispersion should be reduced, resulting in better controlled noise impacts on airport surrounding communities. Current civil aviation guidance systems operate with real time corrective actions to maintain the aircraft trajectory as close as possible to a space-indexed planned trajectory while the flight management system copes indirectly with overfly time constraints. In this paper, we consider the design of new longitudinal guidance laws so that accurate vertical tracking is achieved while overfly time constraints are satisfied. Here, distance to land is adopted as the independent variable for the aircraft flight dynamics since it is today available on board aircraft with acceptable accuracy. A representation of aircraft longitudinal guidance dynamics is developed according to this spatial variable and a space-indexed nonlinear inverse control law is established to make the aircraft follow accurately a vertical profile and a desired airspeed. The desired airspeed is defined by an outer spaceindexed control loop to make the aircraft overfly different waypoints according to a

World air transportation traffic has known a sustained increase over the last decades leading to airspace near saturation in large areas of developed and emerging countries. For example, today up to 27,000 flights cross European airspace every day while the number of passengers is expected to double by 2020. Then safety and environmental considerations urge today for the development of new guidance systems providing improved accuracy for spatial and temporal trajectory tracking.

In the near future air traffic management environment defined by SESAR (Single European Sky ATM Research) and NextGen (Next Generation Air Transportation System) projects, aircraft will progress in four dimensions, sharing accurate airborne predictions with the ground systems, and being able to meet time constraints at specific waypoints with high precision when the traffic density requires it [START_REF] Pappas | A Next Generation Architecture for Air Traffic Management Systemse[END_REF][START_REF] Korn | 4D FMS for Increasing Efficiency of TMA Operations[END_REF]. Then better separation and sequencing of traffic flows will be possible while green climb/descent trajectories will be eased in terminal areas.

Current civil aviation guidance systems generate real time control actions to maintain the aircraft trajectory as close as possible to a planned trajectory provided by the Flight Management System or to comply with ATC tactical demands based either on spatial or temporal considerations [START_REF] Miele | Guidance Strategies for Near-Optimum Takeoff Performance in Wind Shear[END_REF][START_REF] Miele | Optimization and Gamma/Theta Guidance of Flight Trajectories in a Windshear[END_REF]. While wind remains one of the main causes of guidance errors [START_REF] Miele | Optimal Trajectories and Guidance Trajectories for Aircraft Flight Through Windshears[END_REF][START_REF] Psiaki | Analysis of Aircraft Control Strategies for Microburst Encounter[END_REF][START_REF] Psiaki | Thrust Laws for Microburst Wind Shear Penetration[END_REF], these new requirements for improved ATC are attended with relative efficiency by current airborne guidance systems. These guidance errors are detected by navigation systems whose accuracy has known large improvements in the last decade with the hybridization of inertial units with Global Positioning System information.

Nevertheless, until today vertical guidance needs improved performances [START_REF] Singh | Decoupling in a Class of Nonlinear Systems by State Feedback[END_REF][START_REF] Stengel | Toward Intelligent Flight Control[END_REF] and corresponding errors especially during manoeuvre [START_REF] Sandeep | Optimal Nonlinear Estimation for Aircraft Flight Control in Wind Shear[END_REF] may remain large considering the time-based control laws which are applied by flight guidance systems [START_REF] Psiaki | Control of Flight Through Microburst Wind Shear Using Deterministic Trajectory Optimization[END_REF].

In this paper, we propose a new approach to design vertical guidance control laws so that accurate vertical tracking and overfly time are insured. Instead of using time as the independent variable to represent the guidance dynamics of the aircraft, we adopt distance to land, which can be considered today to be available online with acceptable accuracy and availability. A different representation of aircraft vertical guidance dynamics is developed according to this spatial variable where the overfly time becomes a controlled variable. Then a nonlinear inverse control law based-on follow accurately a vertical profile and a desired airspeed [START_REF] Bouadi | Space-Based Nonlinear Dynamic Inversion Control for Aircraft Continuous Descent Approach[END_REF][START_REF] Bouadi | Aircraft Trajectory Tracking by Nonlinear Spatial Inversion[END_REF][START_REF] Bouadi | Aircraft Time-2D Longitudinal Guidance Based on Spatial Inversion of Flight Dynamics[END_REF]. This leads to regulate on a spaceindexed basis the desired airspeed to meet constraints related to the stall speed and the maximum operating speed while making the aircraft to overfly different waypoints according to a planned time-table .   Simulation experiments with different wind conditions are realized for a transportation aircraft performing a direct descent approach for landing under this new guidance law. The resulting guidance performances are compared with those obtained from an equivalent time-based guidance control law.

II. Aircraft Longitudinal Flight Dynamics

The motion of an approach/descent transportation aircraft along a landing trajectory is referenced with respect to a Runway Reference Frame (RRF) where its origin is located at the runway entrance as shown in Fig. 1. The vertical plane components of the inertial speed are such as:

ẋ = -V a cos γ a + w x (1a) ż = V a sin γ a + w z (1b)
where x and z are the vertical plane coordinates of the aircraft centre of gravity in the RRF, V a is the airspeed modulus, γ a is the flight path angle, w x and w z are the wind components in the RRF.

The airspeed and the flight path angle are given by:

V a = ( ẋ -w x ) 2 + ( ż -w z ) 2 (2a) γ a = -arctan ż -w z ẋ -w x (2b)
Adopting classical assumptions such as the RRF being an inertial frame, local flatness of the Earth, constant aircraft mass, the translational acceleration equations can be written as:

mẍ = -T cos θ + D(z, V a , α) cos γ a + L(z, V a , α) sin γ a (3a) mz = T sin θ -D(z, V a , α) sin γ a -mg + L(z, V a , α) cos γ a (3b)
with T , D and L are respectively the thrust, drag and lift forces, α denotes the angle of attack:

α = θ -γ a (4)
The lift and drag forces are given by:

L = 1 2 ρ(z)V 2 a SC L with C L = C L0 + C Lα α (5a) D = 1 2 ρ(z)V 2 a SC D with C D = C 0 + C 1 α + C 2 α 2 (5b)
where ρ(z), S, C L and C D represent the air density with respect to the altitude, the wing surface area, the lift and drag coefficients, respectively. C L0 , C Lα , C 0 , C 1 and C 2 are aerodynamic coefficients. Assuming first order dynamics with time constant τ for the thrust delivered by the engines, we get:

Ṫ = 1 τ (δ th -T ) (6)
Under the above assumptions, the pitch rate dynamics is described by the following state equations:

θ = q with q = M I y (7)
where M , I y denote respectively the pitch moment and inertia moment according to the aircraft lateral axis where M is given by:

M = 1 2 ρ(z)V 2 a ScC m with C m = C m0 + C mα α + C mq qc 2V a + C m δe δ e (8)
where C m0 , C mα , C mq and C m δe are non-dimensional stability derivatives. From Eq. (3a) and Eq.

(3b), state equations of airspeed and flight path angle with respect to airspeed can be rewritten such as:

Va = 1 m T cos α -D(z, V a , α) -mg sin γ a + m ẇx cos γ a -ẇz sin γ a (9a) γa = 1 mV a T sin α + L(z, V a , α) -mg cos γ a -m ẇx sin γ a + ẇz cos γ a (9b)
The time-indexed longitudinal flight dynamics are given by Eq. (1a), Eq. (1b), Eq. ( 6), Eq. ( 7), Eq. (9a) and Eq. (9b). The independent control inputs to the above flight dynamics are chosen here to be the pitch rate q and the throttle setting δ th while w x and w z are perturbation inputs.

B. Space-Indexed Longitudinal Flight Dynamics

Considering that during an approach/descent manoeuver without holdings the distance-to-land time function x(t) is invertible, it is possible to express during these manoeuvers all the flight variables with respect to x and its derivatives instead of time t. It is possible to write for any variable ( * ):

d * dx = d * dt . dt dx = 1 V G . d * dt ( 10
)
where the ground speed at position x and time t is given by:

V G = dx dt = -V a cos γ a + w x (11)
The first, second and third time-derivative of variable (x) are written ( ẋ, ẍ, ... x ), respectively.

Whereas, the k th derivative of variable ( * ) with respect to x is written ( * [k] ) and Eq. (1a), Eq.

(1b), Eq. ( 6), Eq. ( 7), Eq. (9a) and Eq. (9b) are rewritten as:

z [1] = dz dx = dz dt . dt dx = V a sin γ a + w z V G (12a)
θ [1] = q V G (12b)

q [1] = M I y V G (12c) T [1] = δ th -T τ V G (12d) V [1] a = 1 mV G T cos α -D(z, V a , α) -mg sin γ a + m ẇx cos γ a -ẇz sin γ a (12e) γ [1] a = 1 mV a V G T sin α + L(z, V a , α) -mg cos γ a -m ẇx sin γ a + ẇz cos γ a (12f)
with the additioanl equation:

t [1] = 1 V G (13)

III. Vertical Trajectory Tracking Control Objectives

Two main guidance objectives have been considered in this study:

1. To follow accurately a space-referenced vertical profile z d (x) which possibly has been generated by the flight management system, in accordance with economic and environmental objectives and constraints, 2. To respect a desired time table t d (x) for the progress of the aircraft towards the runway or any other reference point.

The desired time table as well as the desired vertical profile should be generated for the whole flight by the flight management system in charge of finding a trajectory meeting satisfactorily economic objectives (expressed in general in terms of fuel consumption and flight time) and air traffic management constraints while flight domain (height and speed limits) constraints must be satisfied.

Trying to meet directly the second objective in presence of wind can lead to hazardous situations with respect to airspeed limits. So this second objective is taken into account through the on-line definition of a desired airspeed to be followed by the aircraft. In this study, it is supposed that accurate online estimates of wind parameters are available [START_REF] Miele | Optimal Trajectories and Guidance Trajectories for Aircraft Flight Through Windshears[END_REF].

From the desired time table t d (x), we get a desired ground speed V G d (x):

V G d (x) = 1/ dt d dx (x) (14)
then, taking into account an estimate of the longitudinal component of wind speed, a spacereferenced desired airspeed V a d (x) can be defined with the consideration of lower and upper limits for the airspeed:

• For lower speeds, the stall speed V S plus a margin ∆V min must be introduced:

V a d (x) = max V S (z d (x)) + ∆V min , V G d (x) -ŵx (x) ( 15
)
where ŵx is the estimate of the horizontal wind speed.

• For higher speeds, the maximum operating speed V M O must be introduced:

V a d (x) = min V M O (z d (x)), V G d (x) -ŵx (x) (16)
• In all other cases:

V a d (x) = V G d (x) -ŵx (x) (17)
Then considering that q and δ th are the control inputs while z and V a are the controlled outputs, we can derive Eq. (12a) and Eq. (12e) with respect to space until we get:

V [2] a = 1 V 2 G A V (z, α, V a , T, W ) + B Vq (z, α, V a , T, W )q + B V T (z, α, V a , T, W )δ th (18a) z [3] = 1 V 2 G A z (z, α, V a , T, W ) + B zq (z, α, V a , T, W )q + B z T (z, α, V a , T, W )δ th (18b)
where W represents the parameters w x , w z , ẇx , ẇz and ẅx , ẅz which can be expressed successively.

The exact expressions of components A V , B Vq , B V T and A z , B zq , B z T in Eq. (18a) and Eq. (18b) are shown in Appendix. 1.

Since the B i terms are in general different from zero, the relative degrees [START_REF] Bouadi | Aircraft Trajectory Tracking by Nonlinear Spatial Inversion[END_REF] of the output variables V a and z are equal respectively to 1 and 2. According to these values, there are no internal dynamics to worry about in the present case since stabilizing the output variables along a smooth trajectory will imply the stabilization of all other state variables.

IV. Space-Indexed Against Time-Indexed Reference Trajectories

In the literature, countless control techniques have been designed for aircraft trajectory tracking using time as the independent variable [START_REF] Magni | Robust Flight Control, A Design Challenge[END_REF] while quite nothing has been published until recently with space as the independent variable [START_REF] Bouadi | Aircraft Trajectory Tracking by Nonlinear Spatial Inversion[END_REF][START_REF] Bouadi | Aircraft Time-2D Longitudinal Guidance Based on Spatial Inversion of Flight Dynamics[END_REF]. However, many ATC solicitations to aircraft guidance can be considered to introduce space based constraints (time separation at a given waypoint, continuous descent approaches, time metered approaches for optimal use of runways, etc). The use of classical time-indexed guidance systems in these situations appears to contribute to the Flight Technical Error (FTE) of the guidance system. Then, to display the interest for this new approach, in this section it is shown how for general aircraft operations linear decoupled space and time-indexed guidance dynamics are not always equivalent.

Considering the relative degrees of and, nonlinear inverse control techniques [START_REF] Bouadi | Aircraft Trajectory Tracking by Nonlinear Spatial Inversion[END_REF][START_REF] Bouadi | Aircraft Time-2D Longitudinal Guidance Based on Spatial Inversion of Flight Dynamics[END_REF] can be used to make these guidance variables satisfy decoupled linear spatial dynamics such as:

2 k=0 a V k V a -V a d [k] = 0 (19a) 3 k=0 a z k z -z d [k] = 0 (19b)
where the corresponding characteristic polynomials According to derivation rules for composed functions, we get:

ξ [1] z = ξz V G (20a) ξ [2] z = 1 V 2 G ξz - ξz VG V G (20b) ξ [3] z = 1 V 3 G ... ξ z -3 ξz VG V G + ξz 3 V 2 G V 2 G - VG V G (20c)
and

ξ [1] Va = ξVa V G (21a) ξ [2] Va = 1 V 2 G ξVa - ξVa VG V G (21b)
with ξ z (x) and ξ Va (x) being the tracking errors related to the desired altitude z d (x) and to the desired airspeed profile V a d (x), respectively:

ξ z (x) = z(x) -z d (x) (22a) ξ Va (x) = V a (x) -V a d (x) (22b)
Then Eq. (19a) and Eq. (19b) can be rewritten as follows:

ξ [2] Va (x) + k 1v ξ [1] Va (x) + k 2v ξ Va (x) = 0 (23a) ξ [3] z (x) + k 1z ξ [2] z (x) + k 2z ξ [1] z (x) + k 3z ξ z (x) = 0 (23b)
where k 1v , k 2v , k 1z , k 2z and k 3z are real parameters such as the roots of s 2 + k 1v s + k 2v and s 3 + k 1z s 2 + k 2z s + k 3z produce adequate tracking error dynamics (convergence without oscillation in accordance with a given space segment) with s denoting Laplace variable.

It appears that when replacing in Eq. (23a) and Eq. (23b) the space derivatives of the outputs by the expressions given by Eq. (20a) to Eq. (21b), we get nonlinear coupled time dynamics for the altitude and the airspeed errors. Only in the case of a constant ground speed where the space and temporal derivatives are related by:

ξ [k] z = ξ (k) z V k G (24a) ξ [k] Va = ξ (k) Va V k G (24b)
we get equivalent linear decoupled time dynamics given by: ...

ξ z + k 1z V G ξz + k 2z V 2 G ξz + k 3z V 3 G ξ z = 0 (25a) ξVa + k 1v V G ξVa + k 2v V 2 G ξ Va = 0 (25b)
This case corresponds to a no wind situation where the airspeed of the aircraft is maintained constant.

In the case where VG remains constant (i.e. V G (t) = V G (t 0 ) + VG .(t -t 0 )) over a time (space) span starting at t 0 , Eq. (25a) and Eq. (25b) become:

...

ξ z + k 1z V G -3 VG V G ξz + k 2z V 2 G -k 1z VG + 3 V 2 G V 2 G ξz + k 3z V 3 G ξ z = 0 (26a) ξVa + k 1v V G - VG V G ξVa + k 2v V 2 G ξ Va = 0 (26b)
then the above decoupled dynamics have time variant parameters and the predictivity (time of response) of these dynamics is lost. It can be however shown that if VG is very small with respect to V G , these dynamics remain stable.

Then the adoption of time based reference trajectories are of interest when guidance requirements can be better expressed with respect to space (especially when time constraints at specific waypoints are considered). Then it appears that adopting in this case a space based trajectory tracking technique should avoid this source of error.

V. Space-Indexed Nonlinear Dynamic Inversion Control Laws

In this section the space-indexed nonlinear dynamic inversion control technique introduced in [START_REF] Bouadi | Space-Based Nonlinear Dynamic Inversion Control for Aircraft Continuous Descent Approach[END_REF][START_REF] Bouadi | Aircraft Trajectory Tracking by Nonlinear Spatial Inversion[END_REF][START_REF] Bouadi | Aircraft Time-2D Longitudinal Guidance Based on Spatial Inversion of Flight Dynamics[END_REF] to perform aircraft trajectory tracking is displayed. The desired vertical profile z d (x) is supposed to be a smooth function of x (in the considered application x is the distance to touchdown) while considering expressions of airspeed constraints in Eq. ( 15), Eq. ( 16) and Eq. ( 17), V a d is supposed to be a piecewise smooth function of x.

Since in general flight conditions the control matrix given by:

Γ =     B zq B z T B Vq B V T     (27)
is invertible [START_REF] Bouadi | Aircraft Trajectory Tracking by Nonlinear Spatial Inversion[END_REF][START_REF] Bouadi | Aircraft Time-2D Longitudinal Guidance Based on Spatial Inversion of Flight Dynamics[END_REF], it is possible to perform the dynamic inversion to get effective trajectory tracking control laws [START_REF] Sandeep | Optimal Nonlinear Estimation for Aircraft Flight Control in Wind Shear[END_REF][START_REF] Psiaki | Control of Flight Through Microburst Wind Shear Using Deterministic Trajectory Optimization[END_REF]. So we get:

    q d (x) δ th (x)     = Γ -1     V 2 G D z (x) -A z V 2 G D Va (x) -A V     (28)
with:

D z (x) = z [3] d (x) + k 1z ξ [2] z (x) + k 2z ξ [1] z (x) + k 3z ξ z (x) (29a)
D Va (x) = V [2] a d (x) + k 1v ξ [1] Va (x) + k 2v ξ Va (x) (29b)
Observe here that while the successive spatial derivatives of desired outputs z d (x) and V a d (x) can be directly computed, the successive spatial derivatives of actual outputs z(x) and V a (x) in Eq.

(29a) and Eq. (29b) can be computed from Eq. (12a) to Eq. (12f) where the wind parameters must be replaced by their estimates. This yields, the control scheme displayed in Fig. 2 where the autopilot operates in the temporal frame to control aircraft attitude while the auto-guidance operates in the space frame.

In order to make the aircraft overfly different waypoints according to a planned time-table

t d (x)
, a simple outer-loop PID controller is introduced. Desired airspeed is computed and regulated to meet constraints related basically to the desired ground speed V G d (x), the minimum allowable speed and the maximum operating speed. Desired ground speed is defined based on the reference time-table t d (x) according to Eq. ( 14). Then the PID speed versus space controller is expressed as:

u t (x) = K p e t (x) + K d de t dx (x) + K i x x-∆x e t (Θ)dΘ ( 30
)
where ∆x is small space interval and:

e t (x) = t(x) -t d (x) (31)
and

V C a d (x) = V a d (x) + u t (x) (32)
The corrected desired airspeed V C a d must also satisfy the airspeed constraints defined in Eq. ( 15), Eq. ( 16) and Eq. ( 17). One of the objectives to be met when applying a space-indexed NDI control to the longitudinal guidance dynamics of an aircraft is to make the vertical tracking error follows dynamics as defined by Eq. (23a) and Eq. (23b) (convergence without oscillation over a given space segment). This objective traduced in the time frame leads to the tracking error dynamics of Eq. (25a) and Eq. ( 25b) which are adopted here to apply the time-indexed NDI control and then perform a fair comparison.

The proposed guidance approach is illustrated using the Research Civil Aircraft Model (RCAM) which has the characteristics of a wide body transportation aircraft [START_REF] Magni | Robust Flight Control, A Design Challenge[END_REF] with a maximum allowable landing mass of about 125 tons with a nominal landing speed of 68m/s. There, the control signals are submitted to rate limits and saturations as follows:

- Here the performances obtained for vertical guidance from the space-indexed NDI and those obtained from the time-indexed NDI are compared in no wind situation. initially late according to the planned time table. It appears clearly that the aircraft increases its airspeed to the maximum operating speed during 12000m until it catches up its delay.

Since except at initial transients the performances look similar, Fig. 6 In the no wind situation, the performances obtained with respect to the desired time table and the airspeed management for a delayed or an early aircraft are displayed. maintains its airspeed at the maximum operating speed until it compensates the initial delay. In the second situation the aircraft is initially in advance with respect to the planned time table and in this case the speed controller sets its airspeed to the minimum allowable speed until the time tracking error is eliminated.

B. Simulation Results in the Presence of Wind

Since in this study the problem of the online estimation of the wind components has not been tackled, it has been supposed merely that the wind estimator will be similar to a first order filter with a time constant equal to 0,35s in one case (time NDI guidance) and with a space constant equal to 28m in the other case (space NDI guidance). Then the filtered values of these wind components have been fed to the respective NDI guidance control laws. Here a tailwind with a mean value of

VII. Conclusion

In this paper a new longitudinal guidance scheme for transportation aircraft has been proposed.

The main objective here has been to improve the tracking accuracy performance of the guidance along a desired longitudinal trajectory referenced in a spatial frame. This objective covers both vertical guidance and compliance with an overfly time table. This has led to the development of a new representation of longitudinal flight dynamics where the independent variable is ground distance to a reference point.

The nonlinear dynamic inversion control technique has been adopted in this context so that tracking errors are made to follow independent and asymptotically stable spatial dynamics around the desired trajectories. It has been shown that a similar tracking objective expressed in the time frame cannot be equivalent when the desired airspeed changes as it is generally the case along climb and approach trajectories for landing. with Υ(z, α, V a , T, W ) and F (z, α, V a , T, W ) are such as:

Υ = 1 V 2 G -2V a Va γa + 2 Va γa (w x cos γ a -w z sin γ a ) -V a γ2 a (
and

B zq = 1 V 2 G 1 m (w x sin γ a + w z cos γ a ) -T sin α - 1 2 ρ(z)V 2 a SC 1 -ρ(z)V 2 a SC 2 α + 1 mV a -V 2 a + V a (w x cos γ a -w z sin γ a ) T cos α + 1 2 ρ(z)V 2 a SC Lα (38a)
B z T = 1 V 2
G cos α mτ (w x sin γ a + w z cos γ a ) + sin α mV a τ -V 2 a + V a (w x cos γ a -w z sin γ a )

In the above equations the temporal derivatives u and ü with u ∈ {x, z, γ a , V a , w x , w z } are related with the spatial derivatives of u by: u = u [1] V G (39a)

ü = u [2] V G 2 + u [1] V [1] G V G (39b)
The deterministic horizontal wind speed component is expressed as:

W x (z) = W 0 (z) ln z z 0 (41a) W 0 (z) = W * 0 cos(ωz + ϕ 0 ) (41b)
where ω and W * 0 denote the circular space frequency and magnitude of the considered wind component.

The stochastic wind components adopt Dryden spectrum model [START_REF] Campbell | A Spatial Model of Wind Shear and Turbulence for Flight Simulation[END_REF] generated from two normalized white gaussian noise processes through linear filters such as:

H δx (s) = σ x 2L xx V a 1 1 + Lxx Va s (42)
and

H δz (s) = σ z L zz V a 1 + √ 3 Lzz Va s 1 + Lzz Va s 2 (43)
Here L xx and L zz are shape parameters (turbulence lenghts) such as:

• For z ≤ 305m:

L xx = z (0.177 + 0.0027z) 1.2 (44a)

L zz = z (44b)
• For z > 305m:

L xx = L zz = 305m (45)
where σ x and σ z represent standard deviations of independent processes such as:

σ z = 0.1W 20 ( 46
)
and W 20 is the horizontal wind speed at 20ft above ground level.

• For z ≤ 305m: σ x = σ z (0.177 + 0.0027z) 0.4 (47)

• For z > 305m:

σ x = σ z (48)
Time and spatial derivatives of the wind components are then given by:

ẇx = W xx ẋ + W xz ż + W xt (49)
with:

W xx = ∂W x ∂x W xz = ∂W x ∂z W xt = ∂W x ∂t (50)
and

ẇz = W zx ẋ + W zz ż + W zt (51)
with:

W zx = ∂W z ∂x W zz = ∂W z ∂z W zt = ∂W z ∂t (52)

V 3 S

 3 a = airspeed, m/s V G = ground speed, m/s γ a = flight path angle w.r.t airspeed, rad θ = pitch angle, rad α = angle of attack, rad L = lift force, N D = drag force, N T = thrust force, N C L = lift force coefficient C D = drag force coefficient m = mass, Kg q = pitch rate, rad/s M = pitch moment, N.m C m = pitch moment coefficient I y = pitch inertia moment, Kg.m 2 ρ = air density, Kg/m = wing surface area, m 2 g = gravity acceleration, m/s 2 τ = engine time constant , s δ th = throttle setting , rad δ e = elevator deflection , rad w x = longitudinal wind component , m/s w z = vertical wind component , m/s

Fig. 1

 1 Fig. 1 Aircraft attitude and acting forces.

2 k=0 a V k s k and 3

 3 k=0 a z k s k are chosen to be asymptotically stable with adequate transients and response times.

Fig. 2

 2 Fig. 2 Proposed flight control structure

Figure. 3 Fig. 3 Fig. 4 Fig. 5 15 - 5 Fig. 6

 33451556 Figure. 3(a) and Fig. 3(b) display respectively altitude tracking performances resulting from time-indexed NDI and space-indexed NDI guidance schemes. While Fig. 4(a) and Fig. 4(b) provide closer views of altitude and tracking performance during initial transients, it appears clearly that in both cases the spatial-indexed NDI trajectory tracking technique provides better results: the spatial span for convergence towards the desired trajectories is shortened by about 2000m while convergence is performed with reduced oscillations. Figures. 5(a) and 5(b) display respectively airspeed tracking performances by space-indexed NDI and time-indexed NDI guidance schemes when the aircraft is

  (a), Fig. 6(b), Fig.7(a) and Fig. 7(b) display respectively the evolution of respectively the angle of attack, the flight path angle, the elevator deflection and the throttle setting during the whole manoeuver. Since the angle of attack remains in a safe domain and the considered longitudinal inputs remain by far unsaturated, this demonstrates the feasibility of the manoeuver.

Figure. 8 Fig. 7

 87 Figure. 8 and Fig. 9 show airspeed and time tracking performances in two cases. The first one considers a delay situation for an aircraft according to a reference time-table where the aircraft

Fig. 8

 8 Fig. 8 Delayed initial situation and recover.

Figure. 11 and

 11 Figure. 11 and Fig. 12 display altitude, airspeed and time tracking performances in the presence of the wind when the actual time table is late and in advance situations according to the reference time table, respectively. It appears that the proposed control technique (space-indexed NDI) keeps its performances shown in the sub-section above.

Fig. 9

 9 Fig. 9 Advanced initial situation and recover.

Fig. 10

 10 Fig. 10 Example of wind components realization.

Fig. 11

 11 Fig. 11 Delayed initial situation and recover with wind.

  w x sin γ a + w z cos γ a ) -V a ( ẅz cos γ a + ẅx sin γ a ) + w x (W zx ẍ + W zz z + Ẇzt ) -w z (W xx ẍ + W xz z + Ẇxt ) -2 V G (-Va cos γ a + V a γa sin γ a + ẇx ) -V 2 a γa -V a ( ẇz cos γ a + ẇx sin γ a ) + V a γa (w x cos γ a + w z sin γ a ) + Va (w x sin γ a -w z cos γ a ) + w x (W zx ẋ + W zz ż + W zt ) + w z (W xx ẋ + W xz ż + W xt ) α -T γa cos α + ρ(z)V a Va SC L -1 2 ρ(z)V 2 a SC Lα γa + mg γa sin γ a -m( ẅx sin γ a + ẇx γa cos γ a + ẅz cos γ a -ẇx γa sin γ a ) -m γa V G -V 2 a γa sin γ a + Va w x -V a (W xx ẋ + W xz ż + W xt )

The tracking performances obtained from the spatial-indexed guidance and the time-indexed NDI guidance have been compared through a simulation study considering a descent manoeuver for a transportation aircraft in wind and no wind conditions. It appears already that the proposed approach should result in improved tracking performances as well as in better trajectory predictability, taking advantage of the enhanced capabilities of modern navigation systems while progress in online wind estimation and prediction are expected.

Appendix

The exact expressions of components A V , B Vq , B V T and A z , B zq , B z T in Eq. (18a) and Eq. (18b) are given by:

and

Appendix

In this study, longitudinal wind is expressed here according to [START_REF] Frost | Wind Shear Terms in the Equations of Aircraft Motion[END_REF][START_REF] Campbell | A Spatial Model of Wind Shear and Turbulence for Flight Simulation[END_REF] as:

where W x (z) and δ x,z (V a , z, t) represent the deterministic and stochastic components of the considered wind, respectively.


Christophe Hurter
email: hurter@recherche.enac.fr

Richard Alligier
email: alligier@cena.fr

David Gianazza
email: gianazza@recherche.enac.fr

Stéphane Puechmorel
email: puechmor@recherche.enac.fr

Gennady Andrienko
email: gennady.andrienko@iais.fraunhofer.de

Natalia Andrienko
email: natalia.andrienko@iais.fraunhofer.de

Wind parameters extraction from aircraft trajectories

Keywords: Wind extraction, least squares approximation, air traffic control, data exploration, visual analytics

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Aircraft fly through the air, and the air flows over the Earth's surface. This simple statement highlights the crucial need to know the winds aloft when one wants to navigate over the Earth's surface in a flying machine. Alternatively, one can also guess how the wind flows simply by observing the trajectories of aircraft relative to the ground. Interestingly enough, air traffic controllers already apply this idea in their everyday work. Experienced controllers can roughly estimate the wind force and direction by observing the aircraft trajectories, and comparing the ground speeds of aircraft flying in different directions: aircraft facing the wind have a lower ground speed than aircraft flying in the opposite direction. This basic idea is at the core of the interactive process proposed in this paper, which allows users to extract the wind direction and magnitude from aircraft radar tracks.

Air Traffic Controllers need accurate wind parameters to perform their activity efficiently. For instance, one can reduce the converging speed of two conflicting aircraft by turning one aircraft so that it will face the wind. The wind impact on aircraft ground speed is also used to slow down or speed up aircraft in order to respect a paced landing sequence, optimizing runway usage (i.e. one landing every 3 minutes). The wind parameters are also necessary to make reliable short/medium term trajectory predictions, so as to avoid trajectory conflicts. With the emergence of new operational concepts [START_REF] Andrienko | Milestone Deliverable D3: The ATM Target Concept[END_REF][START_REF] Swenson | Next generation air transportation system (ngats) air traffic management (atm)-airspace project. Reference Material, External Release Version[END_REF] and automated tools for air traffic management, predicting aircraft trajectories with great accuracy has become more and more critical in recent years. For example, medium-term conflict detection and resolution is very sensitive to trajectory prediction uncertainties [START_REF] Alliot | A statistical analysis of the influence of vertical and ground speed errors on conflict probe[END_REF]. In this context, it is crucial to forecast the wind with accuracy within a prediction window of 15 to 30 minutes, at any point in the airspace. The current meteorological forecasting models do not operate within such timeframes and the best alternative is most probably to use the current wind, assuming that it will remain constant during the time interval of the prediction.

Estimating the current wind numerically still remains a difficult problem, as wind measurements through sensors such as meteorological balloons or radar wind profilers are sparse in both space and time. These wind measurements must be processed by a numerical model and the meteorological wind, pressure, and temperature data is updated every N hours (3 hours, usually).

In this paper, flying aircraft are used as passive wind sensors, with their positions and velocities measured through radar detection. These radar measurements are currently used for air traffic management purposes and are easily available to ground systems, in great quantities. When plotting the ground speed magnitude as a function of ground speed direction, for a number of aircraft, some roughly sinusoidal patterns emerge. These visual patterns are a straightforward result of the wind influence on aircraft trajectories (Figure 1).

In the following, we propose an analytical method to extract the wind from these patterns, using a least squares regression. Instead of focusing on one or two aircraft as in [START_REF] Delahaye | Tas and wind estimation from radar data[END_REF], the idea is to take advantage of the amount of data, and consider categories of aircraft with the same speed characteristics. The wind magnitude and direction is then approximated by applying a least squares method to selected radar tracks, considering all flights in the vicinity of the points where the wind is estimated. A simplified ground speed and wind model is used, neglecting the influence of lateral drift (the drift angle is the angle between the longitudinal axis of an aircraft and its path relative to the ground) on the along-track ground speed. With this approximation, which is justified for aircraft flying at high speeds, the number of unknown variables can be considerably reduced and the model can be linearized.

We also propose an interactive semi-automatic process, where the least squares computation is driven by the user who selects the input data and validates the results through a visual interface.

The remainder of the paper is organized as follows: section 2 details the related works on trajectory exploration and wind extraction, while section 3 introduces the principles of wind extraction from a dataset of aircraft trajectories. The dataset itself is described in section 4. The issues that were raised during some preliminary tests of the automatic extraction method are detailed in section 5. The interactive extraction procedure mixing the least squares algorithm and user filtering and adjustments is presented in section 6. Section 7 shows how the wind dynamics can be extracted and displayed in order to help the user to select time windows when filtering the data. Finally, in section 8, we visually compare our results with Météo-France data. We also discuss the confidence intervals of the least squares method and give some numerical results on the comparison with the Météo-France wind data and the wind computed from downlinked aircraft data. Section 9 concludes the paper and gives some perspectives about further improvements and possible applications of our method.

Related work

Trajectory exploration

The data-flow model in [START_REF] Stuart K Card | Readings in information visualization: using vision to think[END_REF] is widely used to perform data exploration. In this paper, we also use this data flow model to transform raw data (i.e. aircraft records) into visualization with a sequence of transformation steps.

There is a rich bibliography on trajectory analysis in information visualization, in particular on direct manipulation to filter and extract relevant aircraft information [START_REF] Hurter | Fromdady: Spreading aircraft trajectories across views to support iterative queries. Visualization and Computer Graphics[END_REF][START_REF] Hurter | Moleview: An attribute and structure-based semantic lens for large element-based plots[END_REF], density map computation to discover boat trajectory interactions [START_REF] Scheepens | Composite density maps for multivariate trajectories[END_REF], aircraft trajectory schematization [START_REF] Hurter | An automatic generation of schematic maps to display flight routes for air traffic controllers: structure and color optimization[END_REF], trajectory bundling [START_REF] Hurter | Smooth bundling of large streaming and sequence graphs[END_REF], Visual Analytics [START_REF] Andrienko | Visual analytics of movement: An overview of methods, tools and procedures[END_REF], knowledge discovery in databases [START_REF] Giannotti | Mobility, data mining and privacy: Geographic knowledge discovery[END_REF][START_REF] Andrienko | Visual exploration of the spatial distribution of temporal behaviors[END_REF], geocomputations [START_REF] Laube | Progress in movement pattern analysis[END_REF], moving object databases [START_REF] Hartmut | Moving objects databases[END_REF], and detection of landing areas [START_REF] Andrienko | Scalable analysis of movement data for extracting and exploring significant places[END_REF]. However, none of these previous works tried to extract wind parameters from trajectories.

Wind extraction

De Haan and Stoffelen [START_REF] De | Assimilation of high-resolution Mode-S wind and temperature observations in a regional NWP model for nowcasting applications[END_REF] show that high resolution wind and temperature observations can be obtained using on-board measurements made by aircraft and transmitted through the data-link capabilities of the Enhanced Mode-S radars (these data will be detailed in the following). They also show how these measurements can be used as input to meteorological models to improve the short-term and small-scale prediction of wind and temperature. These are significant improvements when compared with the current weather forecasts used in Air Traffic Control (ATC), and they will probably benefit ATC operations in the future.

In the meantime, it is still worthwhile to investigate whether basic radar measurements (position, velocity) are sufficient to estimate the wind. Extracting the wind from radar tracks has already been tried in other works. In [START_REF] Delahaye | Tas and wind estimation from radar data[END_REF], an extended Kalman filter is used to estimate the wind from simulated radar tracks. This method requires two trajectory turns for a single cruising aircraft, or, if two aircraft are considered, one turn per aircraft. The airspeed and the turning rate (in the air) are supposed to be constant for each aircraft. The influence of the wind on the trajectory is modeled by the "wind triangle" stating that the ground speed vector of each aircraft is the sum of its airspeed vector and the wind vector. In [START_REF] Lymperopoulos | Sequential monte carlo methods for multi-aircraft trajectory prediction in air traffic management[END_REF], a multi-aircraft trajectory prediction problem is addressed with sequential Monte Carlo methods, focusing on the inaccuracies related to wind forecast errors. The wind is modeled as the sum of two components: the nominal weather forecast and a stochastic error on this weather forecast. The dimensionality and non-linearites of the multi-aircraft problem lead the authors to introduce a new particle filtering algorithm in order to estimate the error on the wind forecast, discarding standard methods such as Kalman filters. The aircraft airspeeds are assumed to be known. As in [START_REF] Delahaye | Tas and wind estimation from radar data[END_REF], the method is validated on simulated trajectories only. All these approaches are based on specific assumptions about the trajectories (e.g. constant turning rate). As these assumptions cannot be ensured for radar-recorded air traffic data, a different approach is required. Left "wind view", scatterplot with aircraft headings, and their corresponding ground speeds. Sinus shapes emerge which show the wind influence on aircraft ground speed.

Wind extraction principle

Our idea is to take advantage of the amount of data, and to consider categories of aircraft having the same speed characteristics. Contrary to Delahaye et al. [START_REF] Delahaye | Tas and wind estimation from radar data[END_REF], our method is validated with real aircraft trajectories instead of simulated ones, and it relies on more than one or two aircraft. Our system uses Information Visualization (Infovis) techniques [START_REF] Stuart K Card | Readings in information visualization: using vision to think[END_REF] including a scatterplot and the visualization of aircraft ground speed and aircraft heading. With enough data, approximately sinusoidal shapes emerge, one for each aircraft category (or average speed category). Figure 1-right shows recorded aircraft trajectories in 2D with their latitude and longitude. We will refer to this view as the "top view". Figure 1-left shows the same dataset in a scatterplot with the Xaxis showing the aircraft direction relative to the ground, and the Y-axis showing the aircraft ground speed. This view will be referred to as the "wind view". The emerging sinusoidal shape is due to the wind influence on aircraft ground speed. The sine angle shift gives the wind direction (aircraft facing the wind have the lowest ground speed), the amplitude of the sine curve gives the wind speed (the amplitude must be divided by two to retrieve the actual wind speed).

Specific units are used in the Air Traffic Control (ATC) community: aircraft altitudes are given in feet (ft) or Flight Levels (FL). For example, FL350 means 35,000 feet above isobar 1013.25 hPa. Distances are given in nautical miles (NM) with 1N M = 1852m, and speeds in Knots (Kts=NM/h). In order to assess our software with air traffic controllers, we have kept these ATC-specific units in the following, instead of using international units.

Wind influence on aircraft velocities

Figure 1 shows a scatterplot of the aircraft velocities for one day of recorded traffic. The direction followed by the aircraft (ground track angle) is represented on the x-axis and the magnitude of the ground speed on the y-axis. When drawing this graph, we used a 40% transparency setting in order to emphasize the emergence of dense areas: pixels become visible only if many plots are drawn at the same location. As seen in Figure 1, some approximately sinusoidal shapes, stacked one upon the other, emerge from our scatterplot. Apart from the sinusoidal curves, this visualization also shows vertical lines which correspond to the accumulation of pixels representing aircraft with the same direction but different speed. The emerging sinusoidal curves give a visual clue as to how the wind globally impacts the ground speed of flying aircraft. Indeed, for each sinusoidal shape, the lateral distance (along the x-axis) between the x-coordinate of the maximum and minimum ground speed values is 180 degrees (i.e. aircraft facing the wind have the minimum ground speed, and those with the wind behind them have the maximum ground speed). The wind magnitude can therefore be deduced by retrieving the maximum and minimum values of each sinusoidal shape and dividing their difference by two. The wind direction can be directly deduced by considering the direction for which the ground speed is at a maximum (i.e. when the wind is pushing the aircraft).

The fact that several sinusoidal patterns emerge reveals several categories of aircraft. Since it is quite rare that a commercial aircraft flies in circles, one such sinusoid cannot be caused by a single aircraft. Each pattern is due to several aircraft flying in different directions, belonging to a same category regarding their airspeed. Therefore it does make sense to group aircraft which have similar average airspeeds.

In order to highlight the wind influence of aircraft ground speed, we used the multivariate visualization software FromDaDy [START_REF] Hurter | Fromdady: Spreading aircraft trajectories across views to support iterative queries. Visualization and Computer Graphics[END_REF]. Figure 2 represents one day of recorded aircraft trajectories; image 1 shows the aircraft heading on the X axis, and the aircraft ground speed on the Y axis. In order to emphasize the sinus wave on aircraft ground speed, one can filter out low altitude records (which corresponds to aircraft landing or taking off) in image 2, and remove aircraft with a vertical speed (climbing or descending aircraft) in image 3. Low altitude or climbing/descending aircraft do not have a stable airspeed magnitude, therefore the ground speed evolutions cannot fit a sinus shape. These records can be considered as noise, and can be removed. In Figure 3 (the filtered view of aircraft trajectories) the sinus shape shows the wind influence on aircraft ground speed at three different clusters of altitudes (FL100, FL200 and FL300). High speed aircraft fly at high altitude (visible in image 2) and have different sinus shape parameters compared to the two other clusters. Each sinus shape corresponds to specific wind parameters (direction and speed). These images clearly show the wind influence on aircraft trajectories at different altitudes.

Ground speed and wind model

Now that we have seen how the wind characteristics emerge from the visualized data, let us introduce a mathematical formulation of the wind influence on aircraft movement. In this section, we will show that the hypothesis of a sinusoidal curve oscillating around an average speed is not exactly true, and that it requires some simplifications and approximations in the underlying model relating the ground speed to the wind.

Considering a flight i, let us denote V i the ground speed along the track followed by the aircraft, and θ i the track angle counted clockwise from the north reference. T i will denote the true airspeed (TAS), and α i the direction towards which the aircraft is heading, i.e. the angle between the longitudinal axis of the airframe and the North reference. Let us denote W the wind magnitude, and φ the wind direction. In the following, vectors will be denoted in bold font (e.g. W), whereas vector magnitudes will be in normal font (W ).

North

North Ground track The relationship between the true airspeed, the ground speed and the wind is illustrated in Figure 4 and simply expressed as follows, using vector notations:

V i = T i + W (1)
When projecting all vectors on the along-track axis, we obtain:

V i = T i cos(α i -θ i ) + W cos(φ -θ i ) (2)
Our aim is to deduce W and φ from several measurements of (θ i (t), V i (t)) made at different times t in a chosen time interval, using several flights i. In this work, we shall assume that all flights belonging to a same category (i.e. similar performances for the airframe structure and engines) fly at the same average cruising speed. Of course, there will remain some disparities within a same category. Even if all aircraft had the same cruising speed in theory, this hypothesis could only be a statistical one, as airlines may operate their flights differently depending on their cost policies.

So, on one hand, the actual dispersion of speeds within one such aircraft category is expected to be relatively large. On the other hand, for aircraft flying at high speeds, the drift angle α i -θ i is expected to be relatively small: for an aircraft flying at about T i = 450 Kts with a cross-wind of 70 Kts, the error made when considering that cos(α i -θ i ) ≈ 1 is about 1 percent of the aircraft speed.

Taking these considerations into account, the exact model of equation 2 can be simplified by replacing T i cos(α i -θ i ) by a unique speed V , for all flights i belonging to a same category. For a flight i with a velocity V i (t) and a track angle θ i (t) measured at time t by radar detection, we then obtain the following simplified model, considering that W and φ remain constant over the chosen time interval:

V i (t) = V + W cos(φ -θ i (t)) (3)
with V ∈ { Vk |k ∈ {1, . . . , c}} the average speed corresponding to the category of aircraft i.

The average speed V does not have to be known and can be considered as an unknown variable, like W and φ. The only pre-requisite is that every flying aircraft must be assigned to one of the existing classes {1, . . . , c} that must be determined in advance, considering the theoretical cruising speed of each aircraft type. These theoretical cruising speeds are the results of a model ( [START_REF]Base of aircraft data (BADA) aircraft performance modelling report[END_REF]) of the airframe and engine performances provided by the aircraft manufacturers. They are available in the Eurocontrol Base of Aircraft Data (BADA, see [START_REF] Nuic | User manual for base of aircarft data (bada) rev.3.9[END_REF]).

The above model can be linearized by introducing two new variables W X = W cos φ and W Y = W sin φ, and considering that cos(φ-θ i (t)) = cos φ cos θ i (t)+sin φ sin θ i (t):

V i (t) = V + W X cos θ i (t) + W Y sin θ i (t) (4)
As the θ i (t) are numerical values obtained from our radar records, we see that Equation 4 is linear with respect to V , W X , and W Y .

This simplification of the initial model drastically reduces the number of unknown variables. With equation 2, the aircraft headings α i were unknown variables. We had one new unknown variable for each straight trajectory segment of each aircraft. With our approximation for the drift angle, we now have only 3 unknown variables W X , W Y , and V when considering one category of aircraft, or c + 2 unknown variables W X , W Y , V1 , . . . , Vc when considering c categories.

In the following, we see how the least squares approximation method can be applied to extract the wind from an over-determined system of equations resulting from several measurements of (θ i (t), V i (t)) from several flights i, over a chosen time interval.

Least squares approximation of the wind

Let us now consider an airspace volume A over a time interval [t 1 , t 2 ], assuming the wind remains constant within this 4D-volume. Let us consider N measurements {(θ j , V j )|j ∈ {1, . . . , N }} of the track angle θ i (t) and ground speed V i (t)), made at different times t within the interval [t 1 , t 2 ], measured from several flights i.

Let us start with a simple case, and imagine that all aircraft flying through A belong to a single aicraft category of mean speed V . If the quality of the available data is sufficient, that is if we have enough data with a correct dispersion of the θ j values, the unknown variables W , φ, and V can be computed from the N measurements of the track angle and velocity, considering the N corresponding instances of equation 4.

In general, N will be much greater than the number of unknown variables, and our model (eq. 4) will not fit the observed data exactly. Let us introduce ǫ j , the difference between the velocity Vj computed from the model and the observed velocity V j . We now have a system of equations expressing linear relationships between the three unknown variables W X , W Y , and V :

ǫ j = V j -W X cos θ j -W Y sin θ j -V j ∈ {1, . . . , N } (5)
We can use the ordinary least squares method to determine the optimum values of W X , W Y , and V minimizing the quadratic error:

E(W X , W Y , V ) = N j=1 ǫ 2 j = N j=1 (V j -W X cosθ j -W Y sinθ j -V ) 2
This is done by solving the linear system involving the partial derivatives of the error with respect to the unknown variables:

∂E(W X , W Y , V ) ∂W X = 0 ∂E(W X , W Y , V ) ∂W Y = 0 ∂E(W X , W Y , V ) ∂ V = 0 (6)
When the associated matrix is invertible, this system will have solution ( ŴX , ŴY , V ) that minimizes the sum-of-squares error. This solution is meaningful when the matrix is well-conditioned.

Finally the wind is obtained using the following equations, remembering that W X = W cos φ and W Y = W sin φ:

W = ŴX 2 + ŴY 2 φ = arctan 2 ( ŴY ŴX ) (7)

Available datasets

For this study, two datasets of aircraft trajectories were available, the first one containing one day of Mode-C multi-radar records from the Paris area (France), and the second one containing one-day records from the experimental Mode-S radar in Toulouse (South-West of France). The dataset details and the differences between Mode-C and Mode-S radar data are explained in the following section. Meteorological data including wind and temperature on a fixed-size 4D-grid have been collected from Météo-France, for the corresponding days.

Trajectory datasets

Aircraft positions are detected by radars dispatched all over Europe. There are two technologies for aircraft position monitoring: primary and secondary radars. Primary radars use an emitted beam and its corresponding reflection on the aircraft body to compute an azimuth (radar angle) and a distance (echo response time). This kind of radar is passive : no data communication is required between the aircraft and the ground station. Nowadays, the use of primary radars is mostly limited to military applications, where the aircraft are assumed to be potentially non-cooperative.

Secondary radars, widely used in Civil Aviation, emit a beam which embeds a query and then compute an azimuth and a distance thanks to the response beam emitted by the aircraft which embeds specific data. There are different types of secondary radars, depending on the data embedded in the aircraft response. Mode-C data contains the aircraft identity and altitude reports in 100 ft intervals. Elementary Mode-S data contains the aircraft identity, altitude reports in 25 ft intervals, and some basic information (flight status, equipment status). Enhanced Mode-S contains useful additional information such as the aircraft velocities (ground speed, true airspeed, indicated airspeed, Mach number), magnetic heading, roll angle, etc.

All of our radar records were obtained from secondary radars, either from a Mode-C multi-radar system located in Paris (France), or from an experimental Enhanced Mode-S radar located in Toulouse (South-West of France). We used the Paris Mode-C dataset in preliminary experiments (see section 5) to test the least squares approximation method presented in section 3.3. This was facilitated by the presence in this data of some extra information -the aircraft categories relative to airspeed capabilities and operating mode (constant calibrated airspeed, or constant Mach number) -made available by a previous post-processing of this data. These preliminary tests showed some of the drawbacks of the fully automatic wind extraction, and motivated the semi-automatic procedure presented in section 6.1. When the Toulouse Mode-S data became available, we decided to evaluate the interactive wind extraction procedure on this new data, allowing us to compare the results (see section 8) with Météo-France data, and also with a wind computed from on-board measurements of the aircraft ground speed and true airspeed. The two trajectory datasets are detailed below.

Mode-C radar data, Paris area

Our Mode-C records were made over an extended area in the Northern part of France, centered on Paris (see Figure 1). The dataset contains 3,712 trajectories composed of 571,580 points with one point every 15 seconds for each aircraft. The recorded attributes are the location (X, Y, using a Polar Stereographic WGS84 projection centered on 47N 0E), the altitude measured by the difference of pressure with isobar 1,013.25 hPa, the ground speed and ground track angle, and a unique aircraft identifier. This unique identifier is helpful to draw lines in the visualization. Distances are counted in nautical miles (NM), speeds in knots (Kts = NM/hour), and altitudes in feet (ft). The angles are in degrees, counted clockwise relative to the North reference. Since we recorded the data from Paris air traffic control center, the lines representing the aircraft trajectories (Figure 1, right) end at the border of the image which corresponds to the limits of the multi-radar coverage of this center.

Mode-S radar data, Toulouse area

We also investigated a second dataset from a single radar ground station located in Toulouse (South of France, Figure 5). We used this second dataset to validate our wind extraction algorithm since it contains additional Mode-S data. Thanks to the wind triangle principle (see 4), we can compute the wind measured on-board the aircraft. We use this wind computed from downlinked Mode-S data as another source of meteorological data, in addition to the reference wind provided by Météo-France presented in the next section.

Geographically, the dataset covers a circular area of radius 170 NM (315 km) centered on Toulouse. The dataset comprises 1,917 aircraft trajectories, with 169,468 radar reports. The average time span of a trajectory is 25 minutes, with one point every 15 seconds.

Météo-France dataset

The meteorological data provided by Météo-France is a 4D-datagrid (latitude, longitude, isobar altitude, time) containing values of temperature, wind direction and magnitude. Météo-France provided us with two different meteorological datasets corresponding to the recorded days in our two radar datasets.

The 4D-grids are slightly different for these two datasets. For the one corresponding to the Mode-C radar data, the grid at a given time and altitude is composed of 151 rows and 101 columns, and is 587 NM wide and 625 NM high. The altitude varies from isobar 1,013.25 hPa to Flight Level 340 (34,000 feet above isobar 1,013.25 hPa), this range being split into 10 steps. The data is given every 3 hours. The grid location starts from the North east, and continues to the South west of France and covers the whole country. However, in our investigation, we will only consider the area which corresponds to our multi-radar coverage.

The grid corresponding to the Mode-S radar dataset is made of 42 altitude levels, and refreshed every hour. Horizontally, the grid size is 0.1 degrees (in latitude and longitude). The visual and numerical comparisons made in section 8 use this 4D-grid. Note however that this precise 4D-grid is not currently available in the French air traffic control centers: The wind data is still updated every 3 hours. Figure 6 shows the grid of wind magnitude at low and high altitudes. We can observe that at low altitude the wind shows boisterousness whereas the wind gradient is smoother at high altitude. One must be aware that this wind data from Météo-France is actually not the "true" wind experienced by flying aircraft. It is the result of a process of data assimilation and smoothing, using observations from various sensors (sounding balloons, wind profiling radars, etc). Collecting and smoothing this data is a relatively long process (a cycle of several hours), so the resulting wind field may not be up-todate and accurate at the time the user will exploit it. Much can be learned, however, from the comparison of our results with this Météo-France wind model, knowing the amount of scientific and computational effort that is devoted to developing high-quality meteorogical models.

Preliminary tests and algorithm tuning

The preliminary tests presented in this section give us some insights into the behavior of the least squares method, when it is used to extract the wind parameters from aircraft trajectory data automatically. The potential issues are illustrated by several examples. We introduce two quality criteria for the wind estimation. Finally, we summarize the remaining issues with the automatic extraction method, and motivate the choices made when designing the interactive wind extraction procedures presented in section 6.1.

The preliminary tests were run on the Paris Mode-C radar records, because this data was readily clusterized by categories of airspeed and operating mode. This clustering allowed us to assume that all radar records within a same category belonged to aircraft with approximately the same average true airspeed in the cruising phase.

Data quality issues

Basically, wind extraction consists in retrieving the estimated wind at a requested 3D-location P (x, y, z) and at a given time t. This requires the data to be filtered in time and space, so as to select the radar plots in the neighborhood of P (x, y, z, t). The least squares approximation is then applied to the filtered data, as described in section 3.3. Some issues regarding the quality of the input data were encountered during the preliminary tests, when applying this basic automatic procedure. For example, if the filtered aircraft plots are not numerous enough, the automatic extraction can fit a sinus shape, but the wind estimation may not be accurate (Figure 7).

Even if the number of aircraft plots is sufficient, the wind estimation can be incorrect due to a limited angle distribution (Figure 8). If the aircraft record angles have a large distribution, the wind estimation can still be erroneous due to multiple possible solutions (Figure 9).

The automatic wind regression tries to minimize the quadratic error taking into account every aircraft record. This creates an incorrect wind estimation when a cluster of records contains outliers. A few outliers can drastically change the wind estimation parameters (Figure 10).

Data quality criteria, and algorithm tuning

In order to cope with some of the issues presented in the previous section, we designed a quality criterion combining the entropy of the ground speed directions and a threshold on the mean-square error after regression. The entropy criterion ensures that we have a sufficient dispersion of the ground track angles, whereas the meansquare criterion limits the dispersion of velocities around the fitted sinusoidal curve. Good entropy is sufficient to ensure that the associated matrix is well-conditioned.

The entropy is computed as follows, considering the distribution of the ground speed directions among n equal bins partitioning the interval [0 • , 360 • ]:

Ent = - i P i ln P i (8)
with the convention 0 × ln0 = 0, and where P i is the empirical probability (using normalized histograms) that a ground speed direction falls in the i th bin. The entropy value is maximum for uniform distribution and falls when the distribution concentrates. Taking into account these data quality criteria, the wind extraction algorithm is the following:

1. filter the data in time and space, in the vicinity of P (x, y, z, t), 2. assess the quality of the filtered data (entropy criterion), 3. if the entropy criterion is met, apply the ordinary least squares method described in section 3.3 to find a suitable sinus shape, 4. assess the quality of the sinus shape fitting (mean-square error criterion). 5. if the quality criteria (i.e. entropy, and mean-square error) are not met, then return to first step and filter data from a larger neighborhood (up to a given maximum size) of P (x, y, z) and repeat the procedure.

After several tests, we empirically defined acceptable criteria with an entropy value above 1.7 (with 20 bins) and a mean-square error below 0.35.

Summary of our preliminary findings

During our preliminary tests, we found that the introduction of data quality criteria solved the issues illustrated in Figures 78, thus improving the efficiency of our extraction method. However, other issues still remain: presence of outliers, incorrect clustering and multiple solutions. We are convinced that the automatic procedure could be improved again, for example by using robust estimation methods instead of the least squares approximation, or by applying more efficient clustering and filtering techniques to the input data. These improvements are left for future work.

In the current paper, we propose to leave some decisions and choices to the human operator about the assessment of the quality of both the input data and the resulting wind estimates. While experimenting with the automatic procedure, we noticed that inconsistent wind estimations could often and easily be spotted by a human being, by visually comparing the wind estimates in neighboring grid cells. In the next section, we propose an interactive procedure combining the wind extraction algorithm with manual data-filtering and curve adjustments.

Interactive wind extraction

In this section, we describe our process to extract wind direction and speed with an interactive system. In order to extract wind parameters, the user can perform an initial automatic process. He or she defines the number of cells to investigate, and then launches the wind extraction. The software then clusters the aircraft positions into 3D space and time, and tries to fit a sinus shape for each cell. The entropy and error criteria are used to automatically invalidate cells with insufficient angle repartition. The extracted wind parameters are then displayed in small multiples with an arrow in each cell.

Figure 11 shows a top view of the adjusted wind in several grid cells at flight level 350, and also the evolution of wind over time in one of the grid cells. We see in this example that the wind changes orientation from North to West during the day. The computation process lasts 2 seconds with a 4x4 grid and 150,000 records. The software provides geographical grids at 4 different altitudes (FL 250 -FL 450). The arrow indicating the wind direction is displayed with a length proportional to the corresponding wind speed.

The user can then easily recognize inconsistent wind extraction when the direction of arrows does not correspond to the neighboring ones. Wind cannot drastically change direction with neighboring cells. The user can then select each cell and manually adjust extracted wind (i.e. use direct manipulation to adjust a better sinus shape), applying a semi-automatic procedure (see section 6.1). The user can also invalidate the cell if there is not enough data or if a suitable sinus shape cannot be adjusted.

This data validation and adjustment is fast: a few seconds for each cell that needs to be adjusted or invalidated. In order to check if the software and procedure are easy enough to use, we asked an air traffic controller to adjust the extracted wind parameters. Within 2 minutes, he corrected more than 10 cells (the ones with strong wind amplitudes and those with inconsistent wind directions).

Semi-automatic procedure

The interactive procedure allowing the user to adjust the wind in a grid cell is the following:

(a) Filtering stage: filtering is performed in space and time. As explained in section 3.1, aircraft with a vertical speed must be removed. Since aircraft with similar The following three sections describe the views available to the user when adjusting the wind in a grid-cell, how the user interacts with the system and assesses the results of his manual adjustment.

Visualization

Our tool displays two main plots: the "top view" and the "wind view" (Figure 12). The top view displays trajectories with a top visualization: the X-axis shows the longitudes, the Y-axis the latitudes. We also use a color gradient to display aircraft altitude: green shades represent low altitudes and blue shades high altitudes. The top view helps users to observe the selected 2D volume which is used to extract wind parameters. More filtering can be performed with range sliders (Figure 12, lower right part). The wind view displays the same trajectories as shown in the top view but with a different scatterplot configuration. The X-axis shows the aircraft direction (0-360 • ), whereas the Y-axis shows the aircraft ground speed. We display transparent dots that correspond to recorded aircraft parameters. In order to visualize whether dots belong to the same trajectory, we connect them with a line. Since the Y-axis represents aircraft direction, many long horizontal lines appear when aircraft direction changes around 0. In order to remove these visual artifacts and reduce cluttering, we connect trajectory points only if the distance between dots is lower than one third of the scatterplot width.

Interactions

In the interactive procedure described at the beginning of section 6.1, the user must be able to:

• define the extent of spatio-temporal data to be investigated (view filtering),

• adjust a sine curve so that it best fits the filtered data, when the default curve adjusted by the least squares algorithm seems inconsistent,

• assess the wind estimation error and, if necessary, change the spatio-temporal data to be investigated.

Each of these steps is described in detail below.

View filtering

In order to extract wind parameters, the user defines the temporal bounding volume to investigate. We use the same interaction techniques available in [START_REF] Hurter | Fromdady: Spreading aircraft trajectories across views to support iterative queries. Visualization and Computer Graphics[END_REF]. The user leftclicks with the mouse pointer on the top view (Figure 12) to define the center of the selected volume and then manipulates range sliders to define the time range, the altitude range and the latitude and longitude range (Figure 12). When manipulating range sliders, the top and the wind view are automatically updated with the filtered aircraft records. The top view displays the full dataset (to provide data context) but with the selected aircraft shown in color, and the non-selected ones in gray (Figure 12).

User adjustments of the sine curve

Phase adjustment Amplitude adjustment In the next step, the user adjusts the shape of the sine curve shown in gray so that it best fits the visualized aircraft plots in the wind view. The shape of the sine curve is defined by the following formula: f (angle) = Amplitude. sin(angle + Shift)

The user can change its phase (the Shif t angle value) by dragging the sinus shape across the wind view. The user can change the sinus curve Amplitude with the mouse wheel. When the user changes the sine wave parameters, the view updates the corresponding estimated wind speed (Amplitude/2) and direction (Shif t). These parameters are displayed as text values; in addition an oriented arrow shows the wind direction (Figure 13).

User assesment of the estimated wind

In order to assess the validity of the sine wave parameters, the user can display two estimation error metrics (Figure 14). The system computes the distance to the sine curve for every aircraft plot. This distance is displayed with vertical yellow lines which start from the middle of the sine wave and whose length is proportional to the computed distance. In addition, the quadratic error is displayed as a transparent red rectangle.

The user can visually assess if the sine curve fits the aircraft plots correctly by trying to reduce the size of the transparent red rectangle (quadratic error), and by reducing the height of the vertical yellow lines (distance to the sine curve for each aircraft record).

Wind dynamic

In section 6, the illustration of the interactive procedure (Fig. 11) showed the extracted wind at several time intervals. This gave us a succint view of the wind evolution during the day. In addition to this interactive inspection procedure, we propose a more systematic way to investigate the dynamics of the wind. Generally, the goal is to find a compromise between the need to use as many trajectories as possible (to ensure sufficient coverage of the 3D space) and the need to use the most recent data (to avoid using outdated records after the wind has changed). The procedure starts with a user-defined division of the territory into compartments. One possibility is to use a regular grid, either rectangular or hexagonal. A more sophisticated alternative is to use a Voronoi tessellation that reflects the distribution of position records and minimizes trajectory distortions [START_REF] Andrienko | Spatial generalization and aggregation of massive movement data[END_REF].

Next, the user divides the time range of available data into equal time intervals (for example, with a length of 1 hr). For each cell of the territory division and each time interval the automatic wind assessment procedure is applied. In the result, time series of wind speed and direction values are produced. These time series are presented for inspection on time graphs and maps. Figure 15-left shows the spatial distribution of speed dynamics at flight level 350. Similarly, Figure 15-right shows the dynamics of directions. In both maps cell 1x1 is highlighted.

We can observe in Figure 15, that the wind speed was quite stable over the whole time period, with a monotonous change from direction 50 • in the morning to direction 110 • in the evening. According to Figure 15, the wind started to blow in the morning, slowed down at mid-day, and started again in the evening. As a result, it can be recommended to apply the analytical procedures, either fully automatic or interactive, only to the recent data.

If the amount of trajectory data is sufficient, this procedure can be applied separately for different ranges of flight levels. By inspecting outliers on the time graph, it is possible to identify regions that require special attention. In particular, it is necessary to apply manual curve fitting to these cells.

Wind extraction results

In this section, we detail our investigations to validate our wind extraction method. We first compare our results with two other sources of meteorological data. Then, we report informal discussions with air traffic controllers. The results presented here were obtained with the second dataset (Mode-S radar data, South-West of France), which contains additional data that can be used to validate our wind extraction, in addition to the Météo-France data. In order to ease parameter comparisons, we used two different designs to compare wind speeds and directions. Both designs use a color gradient (green to red) to show the error magnitude (Figure 16). In addition, the "speed" design uses a vertical line whose length corresponds to the difference in wind magnitude between our approximation and the Météo-France data. The "angle" design uses two lines: The white line shows the wind direction approximated by our method, whereas the gray line shows the direction of the Météo-France wind.

Visual comparison with Météo-France data

Figure 17 shows the differences between the approximated wind and the reference wind provided by Météo-France, at various altitudes and for three different grid sizes. The three grids at the top of Figure 17 display the automatic wind estimation, without user adjustment. Image 1, below these three grids, shows how the user can adjust the sinus shape to better fit the points, and ignore many outliers. Image 2 shows a cell with a low entropy criteria (not enough data). Image 3 shows a perfect match between the extracted wind parameters and the meteorological data.

We can see that the best results are obtained for the altitudes ranging from flight level 300 to flight level 400. There are much fewer valid cells at very high (above flight level 400) and very low altitudes (below flight level 300), and a few of the remaining valid cells show marked differences with the Météo-France wind. Investigating our dataset, we observed that very high altitudes do not contain enough data to produce accurate wind estimations. Furthermore, at low altitudes, most aircraft are climbing or descending and there are few trajectories at a cruising altitude with a stable airspeed. As the records have been filtered so as to remove climbing or descending aircraft, low altitude cells do not contain enough leveled trajectories, with a sufficient variety of ground speed directions.

Comparing the results for wind magnitude ("speed") and wind direction ("angle") in Figure 17, we can observe that the displayed error (in percentage) is smaller in direction than in magnitude, whatever the altitude.

From this visual comparison, we can conclude that our wind approximation is closest to Météo-France data in the grid cells and at altitudes where there are many aircraft flying in different directions, at their cruising flight level, with a stable airspeed. This was to be expected, considering that the quality of the least squares approximation depends on the quality of the input data.

Numerical comparison with Météo-France and Mode-S winds

Let us now present some numerical results, comparing three different wind values: the wind approximated with the least squares method, the Météo-France wind, and a wind computed from the ground speed and the true airspeed downlinked from the aircraft, and available in the Mode-S data provided by the experimental radar in Toulouse (France). These three different wind values are denoted LS, MTO, and Mode-S, respectively. Figure 18 shows boxplots1 of the differences in wind magnitude -in Knots (left) or in percentage of the reference wind (middle) -and wind direction (right) obtained when comparing the three different wind values (LS, MTO, and Mode-S) pairwise. The results in Figure 18 were obtained considering all flight levels, with a 4x4 grid and a time window of 3 hours.

Considering this pairwise comparison of the three different sources, we can see that they all provide consistent wind values. In Figure 18, the difference between the LS approximated wind and the wind obtained from the two other sources (MTO and Mode-S) is higher than the difference between the wind computed from Mode-S data and the Météo-France wind. This is because Figure 18 shows the results for all flight levels, including the LS wind estimates obtained with low quality input data. Of course, results are better when focusing on flight levels where the input data is of sufficient quality. In practice, the quality of the input data depend on several factors: the altitude range and geographic location, and also the time window chosen when filtering the data. Let us now focus on the upper airspace, at flight levels ranging from 350 to 400, and consider the use of several time windows for the estimation of the wind at a given time (16h, as in Figure 18). Figures 19 and20 show the wind deviations in magnitude, for the chosen altitude range (FL 350 to 400) and for several time windows (30mn, 1h, 2h, 3h). Figure 21 shows the deviations in wind direction. Considering sub-figures 19d, 20d and 21d, where the time window (3 hours) is the same as in Figure 18, we see that the differences in magnitude and direction between the three wind values at flight levels ranging from 350 to 400 are much closer one from the other than when considering all flight levels. number of valid cells (i.e. those satisfying our quality criteria) is smaller for the shortest time windows. This should mitigate the statistical interpretation of the boxplots, and expecially concerning the 1h-time-window, where the apparent bias towards higher values for the deviations of the wind direction might be explained by the small sample size. Actually, a time window of 1 to 2 hours size might be the best compromise if we want the estimates to be sufficiently up-to-date.

The results for altitudes ranging from flight level 300 to 350 are similar to the ones presented here (FL350-400). All these numerical results confirm the conclusion of the visual comparison made in section 8.1. They show the good performances of the proposed wind approximation method when the input data is of sufficient quantity and quality. In practice, the domain of application of our method is the upper airspace, at altitudes ranging from flight level 300 to 400, with relatively dense traffic flying in various directions. At these altitudes, we find commercial aircraft of similar performances, flying at their cruising flight level at approximately the same true airspeed.

Confidence intervals and significance testing

The numerical comparisons in the previous section are relevant only if the confidence intervals associated to the least squares estimations of the wind are of a smaller order of magnitude than the differences with the two other winds (MTO and Mode-S). Assuming we had a very large confidence interval around the wind estimate, and assuming that the the Météo-France and Mode-S wind values fall within this interval, the wind estimates and the differences between the three wind values observed in the previous section would lose their significance.

Let us check the size of the confidence intervals associated with the least squares approximation. In the application described in the paper, a QR decomposition of the design matrix is preferred to the classic normal equations, trading computational efficiency for numerical stability. One benefit of the QR solution is the ease of computation of the terms involved in the expression of the test statistics.

Denoting w = W x , W y , V the vector of unknowns, for each of the components wi , i = 1 . . . 3, of w, a confidence interval at level α is given by:

[ wi -t N -3 α/2 e (Z t Z) -1 ii , wi + t N -3 α/2 e (Z t Z) -1 ii ] (9)
where Z is the design matrix and t N -3 α/2 is the value of the student test statistics for confidence level α/2 and N -3 degrees of freedom. On the data available for the study, a very good adequation between the size of the confidence interval and the quality of the sample assessed using the entropy criterion was observed. Figure 22 shows the size of the confidence interval for each cell of a 4x4 grid, for (W x , W y ) (left) and V (right), in Kts, for altitudes ranging from flight level 350 to 400. We see that the order of magnitude of these confidence intervals is much smaller than the difference between the estimated wind (LS) and the other winds (MTO and Mode-S), which comforts our analysis of the results presented in the previous section.

Air traffic controllers' feedback

In order to assess the potential operational interest of our wind extraction process, we conducted informal discussions with air traffic controllers from Aix ATCC2 . Our goal was to assess the validity of our wind extraction process and to understand how air traffic controllers use wind parameters in their aircraft monitoring tasks.

We performed three interviews. Firstly we asked a set of simple questions:

• "How is the wind important in your daily work activity?"

• "How do you retrieve wind parameters?"

• "How often do you verify wind parameters?"

Secondly, we gave a software demonstration and thirdly we explained our algorithm's rationale.

During the interview, all the controllers confirmed that wind parameters are important to their daily activity, but when we asked them how they retrieved these parameters, their response was not unanimous. They often estimate wind parameters only by looking at the aircraft behavior. "When I compare how aircraft turn when facing north or south, I can assess the wind direction and the wind speed". This estimation is not accurate but sufficient to assess approximate wind parameters. Controllers also have a screen which displays estimated wind parameters provided by Météo-France. This data is displayed in 2 tables with four geographically specific points each: South points (Barcelona, Montpellier, Nice, Ajaccio) and North points (C. Ferrand, Dijon, Lyon, Geneva) with 5 Flight Levels (180-390). These tables are updated every 3 to 6 hours.

The controllers confirmed the validity and interest of our wind parameter extraction method. They also explained that our tool is not designed for air traffic controllers who monitor aircraft but rather for the regulator controller, the one that supervises the traffic regulation and does not have to deal in real-time with aircraft. The regulator controller needs to forecast traffic evolution and therefore our dynamic wind parameter extraction could provide valuable information for this operator.

Conclusion

In this paper, after visualizing the sinusoidal patterns resulting from the wind influence on aircraft ground speeds, we have proposed an analytical method to extract the wind magnitude and direction from the radar tracks of aircraft belonging to various speed categories. A simplified model has been introduced, allowing us to drastically reduce the number of unknown variables and to apply the ordinary least squares method to a linearized problem. The proposed simplification consists in neglecting the effects of lateral drift on the along-track speed for aircraft flying at high speeds. As the performance disparities among flying aircraft (even within a same speed category) are of a greater order of magnitude than the effects of the lateral drift, this approach is justified.

An interactive Visual Analytics system has been developed to demonstrate the results of our automatic approach on recorded radar tracks. Users can explore, validate, or adjust the extracted wind parameters. The wind dynamics can also be extracted from the radar tracks and displayed as time series: Knowing the trends in wind evolution can help the operator in the choice of a time window, when filtering the data before extracting the wind. Filtering the data is a compromise between the quantity of data required to perform the extraction, and its temporal and spatial proximity to the point where the wind is approximated.

The extracted wind has been compared with the Météo-France wind grid, and with the wind computed from Mode-S data (ground speed and true airspeed) downlinked from the aircraft. For this purpose, we used a dataset of radar reports from the experimental Mode-S radar in Toulouse (France). As a result, we have shown that the approximated wind is very close to the wind obtained from the other two sources, at least in airspace volumes where sufficient data is available. We have also discussed the limitations of our method: it requires enough input data, with several aircraft flying in various directions. In addition these aircraft should fly at similar constant airspeeds. In practice, we have shown that this occurs in the upper airspace, where commercial aircraft with similar performances fly at similar cruising flight levels. A good compromise for the time window used to filter the input data seems to be between 1 and 2 hours.

To summarize, our wind extraction method is most efficient when applied to the en-route airspace, at altitudes ranging from flight level 310 to 390, where the cruising commercial traffic is of highest density. Conveniently, such airspace volumes of high traffic density are the ones where air traffic controllers most need accurate wind estimations for their trajectory prediction purposes. Some interviews with air traffic controllers confirmed the interest of our approach, from an operational point of view.

Concerning the perspectives of operational use, one could think of feeding the existing meteorological models with the wind approximated from radar data. In our opinion though, this seems a less promising approach than using the aircraft on-board measurement of the true airspeed and ground speed, downlinked to ground systems via Mode-S. However, collecting and using such data does require fully deployed Mode-S datalink capabilities, and also some additional data processing in order to remove some equipment biases (see [START_REF] De | Assimilation of high-resolution Mode-S wind and temperature observations in a regional NWP model for nowcasting applications[END_REF]). There are good hopes that such wind predictions with enhanced accuracy will be made available in the future, at least in the core traffic areas where Mode-S radars are being deployed. In the meantime, our method could be used as an inexpensive alternative to this approach. It could provide up-to-date wind estimations in dense en-route airspace areas, as a complement to the meteorological wind grid which is currently refreshed every hour at best3 . It could also be useful in geographic areas not covered by accurate meteorological models, or where Mode-S capabilities will not be deployed.

As future work, we plan to try robust estimation methods instead of the ordinary least squares approximation. More extensive numerical experiments could also help to tune our data quality criteria. Another promising path could be to take into account some constraints on the wind field (minimum L 2 norm of the Laplacian), so as to improve our wind estimation. Concerning the interactive procedure, further investigations are in progress, to validate the user performance when adjusting the sinus shape. Finally, we plan to design a specific system to emphasize wind dynamic perception for the air traffic controller.

Figure 1 :

 1 Figure 1: right "Top view", one day of recorded aircraft trajectories over Paris area.Left "wind view", scatterplot with aircraft headings, and their corresponding ground speeds. Sinus shapes emerge which show the wind influence on aircraft ground speed.

Figure 2 :

 2 Figure 2: Wind view with one day of recorded aircraft trajectories. Data filtering to highlight wind influence on aircraft trajectories.

Figure 3 :

 3 Figure 3: The sinus shapes show the wind influence on aircraft ground speed at three different clusters of altitudes.

Figure 4 :

 4 Figure 4: Wind and aircraft velocity. T i , True Airspeed. V i ground speed, θ i track angle, α i aircraft heading angle, W wind magnitude, φ wind direction.

Figure 5 :

 5 Figure 5: Radar data (top view, right) from Toulouse Mode-S experimental radar, and "wind view" (left). In the wind view, records are in transparent black, and those of a same trajectory are connected by a colored line (low altitude records are in green , high altitude in blue).

Figure 6 :

 6 Figure 6: Wind speed data provided by Météo-France at low altitude (below FL10, left), and high altitude (above FL340, right).

Figure 7 :

 7 Figure 7: Not enough data to find a suitable sinus shape.

Figure 8 :

 8 Figure 8: Limited range distribution.

Figure 9 :

 9 Figure 9: Impossibility of finding a wind estimation due to multiple solutions.

Figure 10 :

 10 Figure 10: Cluster of outlier records create incorrect wind estimation.

Figure 11 :

 11 Figure 11: Wind parameter extraction at flight level 350 (top view). The small multiples below show the wind in one cell at incremental times (morning, mid-day, afternoon, night).

Figure 12 :

 12 Figure 12: Interface layout with top view (latitude, longitude), wind view (aircraft speed, aircraft direction) and filters with range sliders used to define the 3D temporal volume to investigate.

Figure 13 :

 13 Figure 13: Using direct manipulation techniques, the user can adjust the sine wave curve location (mouse drag), and amplitude (mouse wheel).

Figure 14 :

 14 Figure 14: After adjustment, the yellow peaks correspond to the difference between the sinus shape and the actual aircraft record. The red rectangles correspond to the quadratic error. The figure on the right shows incorrect sinus shape adjustment, with large errors. The left figure shows a better adjustment, with smaller errors.

Figure 15 :

 15 Figure 15: Wind dynamic extraction with a 4x4 grid at FL 350. Each time series shows the direction or speed evolution over time. Cell 1x1 shows how the wind changed direction and speed during the day: in the morning 50 • , 15 Kts, mid-day 70 • , 5 Kts, in the evening 110 • , 25 Kts.

Figure 16 :

 16 Figure 16: Visual conventions for the display of speed or angle differences.

Figure 17 :

 17 Figure 17: Wind parameter extraction and comparison with Météo-France data, with three different grid sizes. Image 1 illustrates a manual user adjustment. Image 2 shows an invalidated cell. Image 3 shows a valid wind extraction.

  Wind dir. deviations (in deg.)

Figure 18 :

 18 Figure 18: Boxplots of wind deviations -in magnitude (left and middle) and direction (right) -for all flight levels, with a 4x4 grid and a time window of 3h.

Figure 19 :

 19 Figure 19: Boxplots of wind magnitude errors (in Kts), for a 4x4 grid, at flight level 350.

Figure 20 :

 20 Figure 20: Boxplots of relative wind magnitude deviations (in p.c.), for a 4x4 grid, at flight level 350.

Figure 21 :

 21 Figure 21: Boxplots of wind direction errors, for a 4x4 grid, at flight level 350.

Figure 22 :

 22 Figure22: Size of the confidence intervals for the estimated wind (left) and for the estimated average true airspeed, for a 4x4 grid, at flight level 350.

Table 1 :

 1 Number of valid cells, using a 4x4 grid at flight levels 350-400, for time horizons of 30mn, 1h, 2h, and 3h.InFigures 19,[START_REF] Nuic | User manual for base of aircarft data (bada) rev.3.9[END_REF] and 21, we can see how the choice of the time window influences the wind estimates. Clearly, a 30-minute time window is too short. It seems that the best results are obtained with time windows of 2 or 3 hours. However, we must be aware that the boxplots are drawn from data of different sizes. As shown in Table1, the

	Time window size	30mn 1h 2h 3h
	Number of valid cells	5	9	13 12

These boxplots were obtained using the boxplot function of the R environment for statistical computing, with its default settings. The box itself represents the interquartile interval, the bold line is the median, and the whiskers represent either an extremum value or at most 1.5 times the interquartile distance.

ATCC: Air Traffic Control Center.

For the Rapid Refresh (RAP) in the U.S.

Acknowledgements

We would like to thank Christophe Baehr and Météo-France for providing the météorological data used in this study, and the French Air Traffic Services provider (DSNA) for making the radar records available to us. Many thanks also to Serge Roux (ENAC) for collecting and post-processing the Mode-S radar data.


est destinée au dépôt

  Basile Dufay, Sébastien Saez, Christophe Dolabdjian, D. Seddaoui, A. Yelon, D. Menard To cite this version: Basile Dufay, Sébastien Saez, Christophe Dolabdjian, D. Seddaoui, A. Yelon, et al.. Improved GMI Sensors Using Strongly-Coupled Thin Pick-Up Coils. Sensor letters, 2009, 7 (3, Sp. Iss. SI), pp.334-338. 10.1166/sl.2009.1088. hal-00988067


Effect Of Parasitic Capacitance on GMI Magnetic Sensor Performance

Sébastien

  Saez, Basile Dufay, Christophe Dolabdjian, A. Yelon, D. Menard To cite this version: Sébastien Saez, Basile Dufay, Christophe Dolabdjian, A. Yelon, D. Menard. Effect Of Parasitic Capacitance on GMI Magnetic Sensor Performance. Key Engineering Materials, 2013, pp.49(1), 124-127. hal-00988087