Mathieu Roussillo
email: mathieu.roussillo@centralesupelec.fr

Philippe Scouflaire

Sébastien Candel

Benedetta Franzelli

Experimental investigation of soot production in a confined swirled flame operating under perfectly premixed rich conditions

Keywords: Soot, Laser Induced Incandescence, Light Scattering, Rich premixed swirled flame

A detailed understanding of the mechanisms underlying soot production is needed to control particle emissions from practical combustion devices. This requires fundamental investigations of these processes in turbulent combustion systems. This study reports experiments on soot production in a perfectly premixed turbulent flame operating under rich conditions. A swirled ethylene/air flame is stabilized in a confined combustor at atmospheric pressure with full optical access. Quantitative measurements of the soot volume fraction ( f v ) are carried out with the Laser Induced Incandescence (LII) technique, allowing the characterization of the soot spatial distribution and temporal intermittency. While f v is detected by LII mainly in the near-wall region of the burner, light scattering measurements bring complementary information on the presence of particles in the whole chamber. The investigation of a rich perfectly premixed turbulent flame represents an unexplored configuration that provides new insights on soot production in a situation where this process is not dominated by mixing between fuel and air. The data gathered in this situation may be of interest in future developments of modeling methods and their validation for soot prediction.

Introduction

Soot emission is now recognized as an important issue for its adverse effects on climate [1] and human health [2]. The problem is not easy to solve because soot is produced by complex processes and the difficulties are compounded by the fact that most practical devices involve turbulent flames. Numerical simulations are still not capable of precisely predicting soot production [3] and it is important to advance soot modeling to help the control and reduction of soot emission from these systems. In addition to the many studies carried out on sooting laminar flames, it is worth examining soot production in turbulent flows that are widely used in industrial applications. In this respect, much of the technical literature has focused on soot in turbulent jet diffusion flames [4][5][6][7][8]. This case has been extensively investigated with a variety of diagnostics, highlighting the effects of turbulent eddies, strain rate and scalar dissipation on the soot volume fraction f v distribution. Recently, Geigle et al. [9,10] have developed a confined swirled turbulent non-premixed configuration, providing new insights on soot production in a system that has many of the features found in industrial burners. However, all these investigations deal with the non-premixed combustion mode, where soot production is essentially controlled by the quality of mixing between fuel and oxidizer that is mainly governed by the turbulence. The local mixture inhomogeneities dominate the soot production process, complicating the interpretation of the experimental data and numerical modeling of soot production. It was then worth asking whether one could study soot production in a turbulent flame operating under perfectly-premixed rich conditions. This has led to the development of the "EM2Soot" burner, in which a rich premixed swirled turbulent flame is established in an optically accessible combustor [START_REF] Roussillo | Proceedings of ASME Turbo Expo GT2018-76205[END_REF]. The interest of this experimental burner is threefold. First, it allows to independently examine effects of various parameters on soot production such as the equivalence ratio or the flame power, bypassing the complexities linked to the presence of local mixture inhomogeneity [START_REF] Roussillo | Proceedings of ASME Turbo Expo GT2018-76205[END_REF]. Second, this configuration can be used to generate a new experimental database for the validation of soot modeling in turbulent flames. In particular, the perfectly premixed rich flame represents an unexplored situation which complements the experimental database of the International Sooting Flame (ISF) workshop [START_REF]International Sooting Flame (ISF)[END_REF]. To this end, ethylene, the reference fuel for the ISF activities, is here retained due to its high propensity to produce soot and the large availability of models for sooting ethylene flame simulations. Third, the investigated rich premixed flames provide information that may be of interest in the design of combustors employing Rich-Quench-Lean (RQL) concepts to reduce pollutant emissions, in which the first stage combustion process is carried out in a rich regime at an equivalence ratio φ 1.8.

In the present work, soot production in the EM2Soot configuration is characterized and this is, to our knowledge, a first of this kind. The experimental setup is briefly introduced in Section 2, together with the Laser Induced Incandescence (LII) and Laser Light Scattering (LS) techniques employed to investigate soot volume fraction and soot particle presence, respectively. Soot production is then characterized in Section 3, by examining the volume fraction distribution, and the spatially layered structure of soot particles. Using light scattering imaging it is also shown that soot particles are present in some regions where their volume fraction is below the detection threshold of LII. Section 3 then discusses the temporal variability of soot volume fraction induced by the turbulent eddies.

Experimental Setup

A detailed description of the EM2Soot burner may be found in [START_REF] Roussillo | Proceedings of ASME Turbo Expo GT2018-76205[END_REF] and is only summarized in what follows. The burner, schematically represented in Fig. 1a, includes three components: the plenum, the injector and the combustion chamber. The chamber features a square cross section and its dimensions are 118 × 118 × 250 mm 3 . An ethylene/air mixture at ambient temperature prepared in the upstream manifold is conveyed to a plenum feeding an injection channel equipped with a radial swirler (Swirl number S = 0.7). It is formed by four quartz windows (90 mm × 250 mm × 3 mm) mounted between metallic corner pieces. One of the corner bars is equipped with 8 thermocouples measuring the temperature on the inner side of the chamber wall. The burner can be operated in a large domain of equivalence ratios and mass flow rates, but this study focuses on the reference operating conditions summarized in Table 1, chosen for the high soot production of the corresponding flame. In this configuration, soot generation is strongly dependent on the wall temperature [START_REF] Roussillo | Proceedings of ASME Turbo Expo GT2018-76205[END_REF].

It is therefore ensured that all experiments are carried out at a fixed initial wall temperature (T = 573 K measured by one of the thermocouples at a height above the burner HAB=113 mm), which is obtained by preheating the chamber with a lean flame to prevent quartz obscuration by soot deposits. During each experimental run, corresponding to a duration of about 30 s, the wall temperature varies by less than 3% so that a quasi-steady state can be assumed for the thermal environment. In addition, using a short acquisition time the soot deposit on the wall can be considered to be negligible during each experimental run. The broadband emission of light from the reference flame, recorded by a Nikkon D700 camera, is displayed in Fig. 1a. Light is emitted from a large volume of the chamber, but the largest intensities originate from the lower part, where a highly luminous yellow flame is observed. One may then deduce that soot concentration is mainly present in the vicinity of the chamber bottom plane. However, the detected signal also depends on the soot temperature through Planck's law so that the intensity only provides a qualitative indication of the presence of soot.

It is then more instructive to deduce 2-D soot volume fraction fields from LII techique since the LII signal is effectively proportional to soot volume fraction [START_REF] Melton | [END_REF].

The detection performance of this technique depends on the laser characteristics and on the optical setup, allowing to capture the presence of soot only when f v exceeds a critical threshold value (here = 7.5ppb). Light scattering (LS) imaging is also used in the present investigation as a complementary optical diagnostic. Its physical interpretation is less easy since the scattered light signal is not directly proportional to f v , being strongly dependent on particle morphology and sizes [14]. For this reason, the LS images are essentially used as an indication of soot presence and of the spatial layout of soot. As the sensitivity of LS is not linked to the volume fraction value, it can be used to detect the presence of soot particles in regions where f v < . LII measurements are performed using a laser sheet focused on the burner central axis (convergent lens f=+1000, cylindric lens f=-50), obtained from a Nd:YAG laser (1064 nm, 10 Hz). The sheet is 0.35 mm thick and extends over 7 cm in the vertical direction. It should be noticed that a uniform fluence is guaranteed only in the central 4 cm of this laser sheet. The laser fluence is set to approximately 0.45 J/cm 2 with a shot duration of ≈ 9 ns. This value belongs to the "plateau" region [15] of the LII signal for this configuration [START_REF] Roussillo | Proceedings of ASME Turbo Expo GT2018-76205[END_REF]. LS measurements are performed using a 9 cm wide laser sheet generated with a Nd:YAG laser doubled to 532 nm with a shot duration of ≈ 9 ns. A low energy fluence (0.09 J/cm 2 ) is imposed to avoid soot incandescence.

The LII and LS signals are captured with two intensified CCD cameras (Princeton, PI-MAX 3, f=10Hz) and stored into a 3-D matrix:

I(x, y, t) ∈ R N x ×N y ×N t +
, with N x = 1024, N y = 1024 and N t = 300, defining a pixel size of 0.1mm. Close-up views are also obtained with a pixel size of 0.05mm by inserting an extension ring to the camera objective. Both cameras are equipped with Nikkor 50 mm f/2 lenses. A bandpass filter centered at 425 nm (50 nm FWHM) is used for LII to reduce parasitic C 2 Swan band emission and flame luminosity [16], whereas a bandpass filter centered at 532 nm is used for LS. A reliable proportionality between LII signals and soot volume fraction is achieved by suitably choosing the acquisition duration of the LII signal [START_REF] Melton | [END_REF]. The signal is recorded with no delay with respect to the end of the laser pulse over a period of 100 ns to maximize the signal to noise ratio. The LS signal is acquired for 20 ns without delay with respect to the laser pulse. The LII signal has been calibrated via the Modulated Absorption Emission technique (MAE [17]) on a laminar diffusion flame [START_REF] Franzelli | Proc. Combust. Inst[END_REF], fixing the value of the soot complex refractive index E(m, 645 nm) = 0.38. It is known that uncertainties in the value of E(m) may lead to a relative error on the soot volume fraction that may be up to a factor two [START_REF] Franzelli | Proc. Combust. Inst[END_REF].

Characterization of soot production

Soot structure

Images of the LII signal in the form of instantaneous soot volume fraction fields along the chamber are displayed in Fig. 1b. It is worth noting that the laser beam enters the chamber from the left side. A net reduction in the intensity level can be observed on the right side, providing an indication on the soot absorption of the laser sheet as it travels in the chamber. Measurements in the right side of the flame are biased and the left side results are only considered in the following. Examining the LII signal, one notices that large volume fractions of soot are mainly present near the wall 1 Images indicate that the f v values are generally below the LII detection threshold in the flame central region. Instantaneous f v values can reach about 0.1 ppm near the wall (|r| > 3 cm) close to the bottom plane (HAB< 4 cm), which is fairly small compared to the levels classically observed in turbulent non-premixed flames [4,7]. Further downstream, the level of instantaneous f v generally decreases. Soot volume fraction appears in the form of layers of approximately 0.5-2.5 mm thickness, which are clearly distinguishable and comparable to previous 1 The immediate vicinity of the wall is not accessible due to the limitation of the optical access, reducing the field of view by about 5 mm for each side. experiments [4,6]. A close-up view of these ligaments is presented for two HABs in Fig. 1c. The effect of turbulent eddies on the ligaments is quite visible, as they are wrinkled and distorted by the turbulent flow. Close to the injection plane (bottom image in Fig. 1c), larger eddy structures can be observed, so that soot particles occupy a reduced surface area, i.e. a high preferential concentration of soot particles prevails. Downstream (top image in Fig. 1c), large eddy structures are less visible reducing the preferential concentration of soot.

One may note that the spatial soot organization differs from that observed in the swirled non-premixed flame investigated by Geigle et al. [10,[START_REF] Geigle | [END_REF], which to our knowledge represents the only experimental database available in the literature for confined swirled sooting flames. This is most likely due to the different combustion modes. In particular, soot appears near the chamber walls in the present configuration in contrast with the pattern obtained in the DLR flame where f v reaches its largest levels in the middle of the chamber and is nearly negligible in the outer recirculation zone close to the walls. Soot ligaments are also shorter than those found in the present experiment but both configurations feature a high spatial preferential concentration of soot, also designated as spatial intermittency. In the DLR flame, soot presence is related to the mixture fraction field and is linked to locally rich zones promoting soot production. The spatial distribution of soot then depends on turbulent mixing. In the present premixed flame, soot production is not governed by the local equivalence ratio since it is homogeneous in the whole chamber. Here, soot is expected to be mainly affected by the turbulence through fresh and burnt gas mixing, its effects on local flame properties (such as temperature or soot precursor concentration) and the stretch induced on soot ligaments by turbulent eddies.

Spatial soot structures are now investigated by exam-ining LS images [6] to collect information on soot presence and complement the LII data. This is exemplified with quasi-simultaneous measurements of LS and LII performed with the same laser sheet generator with two Nd:YAG lasers doubled to 532 nm. Light scattering is measured with a first low energy laser shot 20 ns before the LII shot to prevent any parasitic LII signal and vaporization of smallest soot particles. Results are displayed in Fig. 2. The LII signal (red) only features a few ligaments close to the wall as it pertains to sufficiently high detectable values of f v . The LS signal (black) is superimposed on the LII signal where this signal is present but also recovers more soot ligaments everywhere in the chamber. These ligaments may be constituted by a small number of large aggregates or a big number of nuclei particles, so that the total soot volume fraction is negligible and therefore not detectable by the LII. At this point, it is not possible to know which is the correct interpretation. However, since light scattering of spherical particles in the Rayleigh regime is proportional to d 6 p , the first proposition, i.e. small number of large aggregates with bigger primary particle diameter d p compared to soot nuclei, seems more likely but there is no clear proof of it at this stage. This point might be resolved in the future by performing time-gated LII measurements [START_REF] Franzelli | Proc. Combust. Inst[END_REF]. A 2-D instantaneous field of light scattering is displayed in Fig. 3a. by combining two uncorrelated images. Long soot ligaments highly wrinkled by the turbulent eddies are detected in the whole region. By comparing LII and LS instantaneous results (Figs. 1b and3a, respectively), it can be deduced that sufficiently high values of f v prevail along the walls, a. but soot particles are actually present everywhere in the chamber.

Since soot particles are characterized by a nearly zero Stokes number, the soot pattern detected via LS contains information on the turbulent velocity field as one may reasonably assume that diffusion and thermophoresis effects are negligible. Therefore, the Particle Imaging Velocimetry (PIV) methodology is here applied to 150 pairs of soot LS fields separated by 75 µs using an adaptative PIV algorithm with a minimal interrogation areas of 32x32 pixels. The images have been obtained using a 2x400 mJ Nd:YAG laser doubled at 532 nm (10 Hz, Continuum, Powerlite SL3-PIV) and a Flow sense Eo (Dantec Dynamics) camera equipped with a 532nm filter (10nm bandwidth) to enhance the signalto-noise ratio. The resultant axial velocity U field presented in Fig. 3b only provides qualitative information on the velocity field since, being subjected to mass and number variations, soot particles cannot be considered as "perfect" tracers. It is however possible to recognize the presence of a wide central inner recirculation zone (IRZ), identified by the black isocontour U = 0. Compared to a classical swirled flame, one notices the absence of outer recirculation zones close to the chamber corners. This atypical flow structure is characteristic of a toroidal swirled flame, that expands in the lateral direction and occupies a region located near the combustor bottom plane [START_REF] Degeneve | Proceedings of ASME Turbo Expo GT2018-76152[END_REF]. Unfortunately, this cannot be confirmed by OH* chemiluminescence since the corresponding signal is dominated by black-body radiation from the sooting flame captured by secondary rebounds of the filter [START_REF] Roussillo | Proceedings of ASME Turbo Expo GT2018-76205[END_REF]. However, from the OH* field obtained in the EM2Soot configuration for the richest non-sooting premixed flame, it is possible to identify a toroidal swirled flame [START_REF] Roussillo | Proceedings of ASME Turbo Expo GT2018-76205[END_REF], confirming results obtained for the axial velocity field and the inferred flame structure. A more precise investigation of these is however out of the scope of this work.

Time-averaged and RMS fields of f v

The time-averaged field of soot volume fraction f v is obtained by summing 300 images. Results are shown in Fig. 4a. As already indicated, the asymmetry of the results is due to the soot absorption of the laser sheet that enters the chamber from the left. Time-averaged results confirm that f v is almost negligible in the middle of the chamber and takes detectable values close to the walls. The time-averaged soot volume fraction initially reaches a maximum value of 10ppb (region A for 1 < HAB < 5 cm), then decreases downstream at around 6ppb (region B for 5 < HAB < 12 cm) and finally increases again at 8ppb (region C for 12 < HAB < 18 cm). a. The trend observed for f v differs from what may be concluded from the flame luminosity image 2 in Fig. 1a. By considering the flame luminosity, one can deduce that soot is present in the whole chamber while LII measurements only detect high soot volume fraction levels along the lateral wall. However, the difference is due to the fact that the light intensity signal is integrated along the line-of-sight and cannot be compared to the planar 2-D field of f v obtained with LII. Second, from Fig. 1a, f v is expected to reach its maximum value near the bottom of the chamber and to rapidly decrease from that point. These conclusions are not confirmed by the LII measurements since f v plotted in Fig. 4a features high values at HAB ≈ 16 cm. This definitely indicates that flame luminosity cannot be used to estimate the spatial evolution of f v in the chamber since it is also dependent on the forth power of the temperature. Results from the standard deviation of f v , σ, are displayed in Fig. 4b. One finds that σ is maximum in region A and decreases in regions B and C, whereas f v has a non monotonic behavior. This indicates a complex behavior of soot production that may be investigated by characterizing the temporal intermittency of the soot distribution.

Temporal soot intermittency and turbulence effects

Soot production in turbulent flames is known to be a highly intermittent phenomenon [6,[START_REF] Qamar | [END_REF]22], which may be characterized by the intermittency index Ω [23]:

Ω(x, y) = 1 - 1 N t N t t=1 I(x, y, t) (1)
2 Due to the exposure time for flame luminosity visualization (2µs), Fig. 1a contains time-averaged information. where I(x, y, t) is the presence index:

I(x, y, t) = 0 if f v (x, y, t) < 1 if f v (x, y, t) ≥ (2)
The 2-D intemittency index, displayed in Fig. 5, reaches its minimum in region C, indicating that soot presence is sporadic for small HAB (regions A and B), whereas soot particles are more likely to be continuously detected at higher HAB (region C). The temporal intermittency can strongly affect the interpretation of soot production trend, especially when it is based on timeaveraged fields. Therefore, a weighted time-average operator • is here introduced to take into account the intermittency of soot particles, by weighting time-averaged soot volume fraction by the presence of soot particles:

f v (x, y) = N t t=1 f v (x, y, t)I(x, y, t) N t t=1 I(x, y, t) (3)
Results for the weighted time-averaged field f v are shown in Fig. 5b. Compared to the f v field, the weighted mean f v is higher close to the bottom plane (region A) and slowly decreases downstream towards a nearly constant value. Evolutions with axial position of f v , Ω, f v and σ variables in the near-wall region are summarized in Fig. 6. Here, each variable has been spatially averaged over a surface of 1 cm × 1 cm (r ∈ [-4.9; -3.9] cm). The temporal and spatial evolution of soot volume fraction can be characterized by combining results from Fig. 6 with local information on the probability to detect a certain value of soot volume fraction. For this, the probability density function (PDF) of f v for three axial positions corresponding to regions A, B and C are shown in Fig. 7. Each histogram represents the distribution of instantaneous soot volume fraction in a square of 5 × 5 mm 2 centered on r = -4.5 cm accumulated over 300 images (HAB=3,10,17 cm). In each histogram, f v and σ values are reported. The red line provides an estimation of the PDF for the f v values below assuming a log-normal distribution. Region A is characterized by a high soot temporal intermittency, indicated by the height of the first vertical bar. When soot is detected, the f v threshold is the most probable value and, at the same time, the highest values of f v (about 0.15 ppm) are observed in this region, leading to a maximum σ value in the center of region A (HAB= 3.5 cm). This region corresponds to the IRZ so that this behavior possibly indicates that the large eddies of the turbulent flow observed in this region in Fig. 1c induce a high variability of instantaneous f v values. Thus, soot production can be enhanced or diminished by their actions, a feature discussed in [24,25] for a single vortex interacting with a sooting laminar flames. Downstream, the region close to the wall for 5 < HAB < 12 cm (zone B) features nearly the same level of intermittency as in region A. Its most probable f v value is close to but σ is reduced so that the final time-averaged and weighted time-averaged values are smaller than in region A. Finally, for higher HAB (region C) soot temporal intermittency is drastically reduced. Since σ is the smallest, f v value is high again at HAB= 16 cm. In summary, rich premixed swirled flames feature a strong relation between the spatial structure (Fig. 1c) and the temporal variability of f v . Close to the chamber bottom corners, long ligaments of soot are wrinkled by turbulence forming large eddy soot structures with a spatial preferential concentration which temporally results in a high level of intermittency. In addition, strain rate and curvature effects induced by vortices are known to decrease or enhance soot production, respectively, increasing the variability of f v value along a soot ligament [25]. Downstream, f v is detectable via LII close to the wall and a quasi-homogeneous spatial distribution of soot is observed (Fig. 1c), probably due to an attenuation of the eddies strength near the wall. This results in a thinner PDF for f v and a higher probability to detect soot as a function of time. One may wonder whether soot oxidation plays a role in this configuration. To answer this question an attempt was made to detect OH radicals with planar laser induced fluorescence (PLIF) with a dye laser tuned to λ = 282.92 nm and detection in the band 308±5 nm (ZBPA310 Asahi bandpass filter). No signal could be detected along the height of the burner, indicating that OH volume fraction was negligible. This is in agreement with the fact that at φ = 2.1 the OH concentration is expected to be quite small compared to that prevailing under stoichiometric conditions. Based on the analysis of the flame structure proposed in Section 3.1, one may consider that most of the OH concentration will be localized near the combustor bottom plane. In the absence of this radical and of oxygen (due to rich premixed con-ditions), soot oxidation is expected to be quite limited in this configuration [START_REF] Roussillo | Proceedings of ASME Turbo Expo GT2018-76205[END_REF].

Conclusion

Soot production in turbulent flames is investigated by making use of turbulent combustion operating in the premixed mode. Experiments are carried out on a novel laboratory scale model combustor that is used to stabilize an ethylene/air rich swirling turbulent flame, established in a confined environment allowing optical access to the process. Production of soot is examined using planar laser induced incandescence (LII), providing quantitative measurements of soot volume fraction f v along the height of the burner and planar laser light scattering (LS) by soot particles providing information on the presence and spatial geometry of soot layers. The LII reveals two important trends: soot is mainly present along the wall of the burner highlighting the role of the confinement on soot production in rich premixed flames. The evolution of the mean soot volume fraction along the height of the burner exhibits two locals maxima. The first maximum is located in the lower part of the chamber while the second is formed in the second half of the combustor at a larger distance from the injection plane. A statistical analysis of the two regions reveals that these two maximum values conceal two different trends regarding soot intermittency. Soot intermittency is highest in the lower part of the chamber and decreases along the height of the burner. A weighted mean is proposed to take into account the intermittency of soot particles leading to a rather monotonic evolution of this weighted mean along the height of the burner. Laser light scattering images combined with LII images indicate that soot particles are also present in the central region of the burner but that f v reaches its highest values in the near wall region of the burner. The volume fraction levels in the central region are outside the dynamic range of the LII and their values are too low to be detected.

Figure 1 :

 1 Figure 1: (a) Schematic of the EM2Soot burner. (b) 2-D mapping of instantaneous soot volume fraction by LII as a selective collage of images recorded at different heights (pixel size ∆x = 0.1 mm). (c) Close-up views of soot structure at two HABs: the camera is equipped with an extension ring (∆x = 0.05 mm).

Figure 2 :

 2 Figure 2: Simultaneous measurements of soot presence by LII (red) and LS (black) at λ = 532 nm.

Figure 3 :

 3 Figure 3: a) Image of light scattering intensity, indicating the instantaneous soot presence (combination of two images recorded at different instants). Parasitic background signals from light reflection on the chamber walls are not eliminated as they provide information on the chamber geometry and on the position of the bluff-body (grayscale). b) Axial velocity field derived from PIV treatment of soot light scattering fields. The black contour corresponds to null axial velocity and indicates the IRZ position.

Figure 4 :

 4 Figure 4: 2-D mapping of (a) time-averaged and (b) standard deviation of soot volume fraction from LII.

Figure 5 :

 5 Figure 5: 2-D mapping of (a) temporal intermittency and (b) weighted time-averaged soot volume fraction field from LII.

Figure 7 :Figure 6 :

 76 Figure 7: Evolution of the distribution of the soot volume fraction in the three zones A, B and C with the corresponding log-normal fitting curve.

Table 1 :

 1 Reference flame operating conditions: equivalence ratio, flame power, fuel and air mass flow rates, Reynolds number.

	φ	P tot	ṁF	ṁair	Re
		(kW) (g/s)	(g/s)
	2.1 15	0.319 2.23 20000
	2

Acknowledgments

Support from G.Legros and J.Bonnety (UPMC) for the LII calibration through the MAE technique is gratefully acknowledged. This study is supported by the Air Liquide, CentraleSupelec and CNRS Chair on oxycombustion and heat transfer for energy and environment and by the OXYTEC project, grant ANR-12-CHIN-0001 of the French Agence Nationale de la Recherche.


Zhengyi Wang
email: zhengyi.wang@alumni.enac.fr

Man Liang

Daniel Delahaye

A Hybrid Machine Learning Model for Short-term Estimated Time of Arrival Prediction in Terminal Manoeuvring Area

Keywords: Air Traffic Management, 4D Trajectory Prediction, Clustering, Multi-Cells Neural Network, Machine Learning, Data mining

4D trajectory prediction is the core element of future air transportation system, which is intended to improve the operational ability and the predictability of air traffic. In this paper, we introduce a novel hybrid model to address the short-term trajectory prediction problem in Terminal Manoeuvring Area (TMA) by application of machine learning methods. The proposed model consists of two parts: clustering-based preprocessing and Multi-Cells Neural Network (MCNN)-based prediction. Firstly, in the preprocessing part, after data cleaning, filtering and data re-sampling, we applied principal Component Analysis (PCA) to reduce the dimension of trajectory vector variable. Then, the trajectories are clustered into several patterns by clustering algorithm. Using nested cross validation, MCNN model is trained to find out the appropriate prediction model of Estimated Time of Arrival (ETA) for each individual cluster cell. Finally, the predicted ETA for each new flight is generated in different cluster cells classified by decision trees. To assess the performance of MCNN model, the Multiple Linear Regression (MLR) model is proposed as the comparison learning model, and K-means++ and DBSCAN are proposed as two comparison clustering models in preprocessing part. With real 4D trajectory data in Beijing TMA, experimental results demonstrate that our proposed model MCNN with DBSCAN in preprocessing is the most effective and robust hybrid machine learning model, both in trajectory clustering and short-term 4D trajectory prediction. In addition, it can make an accurate trajectory prediction in terms of Mean Absolute Error (MAE) and Root Mean Squared Error (RMSE) with regards to comparison models.

Introduction

4D trajectory prediction is the core element of future air transportation system [START_REF] Faa | Common trajectory prediction structure and terminology in support of SESAR & NextGen[END_REF]. During the development of Trajectory Based Operation (TBO) concepts of Single European Sky ATM Research (SESAR) and Next Generation Air Transportation System (NextGen) programs, to improve the predictability of air traffic is a crucial and fundamental research topic. 4D trajectory prediction refers to the calculation and prediction of longitude, latitude, altitude and time on the future waypoint sequence based on the historical flight data. According to the time scale, 4D trajectory prediction can be divided into two categories [START_REF] Guan | A strategic flight conflict avoidance approach based on a memetic algorithm[END_REF]:

1. Tactical (short-term) trajectory prediction: A prediction in a short period within several minutes or even shorter. Since the prediction scale is relatively small, minor change may have significant impact on prediction results. Therefore, tactical trajectory prediction requires as much information as possible.

Flight-related information contained in radar or ADS-B data is usually taken; 2. Strategical (long-term) trajectory prediction: A kind of prediction before departure based on the flight plan, which provides the prediction from a macroscopic view. It is mostly applied to fuel consumption and airspace flow evaluation.

4D trajectory prediction can be influenced by numerous factors, such as aircraft weight, pilot actions, wind and temperature. These uncertainties will not only make it difficult to improve the prediction accuracy, but also will decrease the prediction process efficiency as the prediction time becomes longer [START_REF] Tastambekov | Aircraft trajectory forecasting using local functional regression in sobolev space[END_REF]. In this paper, we propose a novel short-term trajectory prediction model, which combines different machine learning techniques, to address the problem of 4D trajectory prediction in Terminal Manoeuvring Area (TMA). The model is mainly divided into 2 parts: preprocessing model and prediction model. The preprocessing model extends and improves the trajectory clustering method proposed by [START_REF] Gariel | Trajectory clustering and an application to airspace monitoring[END_REF], plays a role of preparing data with high quality for the prediction part. In addition, it supports an efficient and accurate prediction process. It contains the following steps: data cleaning, filtering, re-sampling, Principal Component Analysis (PCA), clustering and model training. In the prediction model, a novel technique named Multi-Cells Neural Networks (MCNN) is proposed. It can make parallel trajectory prediction for various shapes of traffic patterns. With Neural Network (NN) for each prediction cell, an accurate prediction for the whole complex trajectories in TMA is realised. Thus, the main contributions of this paper are threefold:

1. A novel hybrid 4D trajectory prediction model combining clustering-based preprocessing and MCNN is developed. A general guideline for selecting the model structures is presented in detail. 2. The proposed model is robust with regards to the data source. The 4D trajectory in a specified airspace containing noise and abnormal trajectories can be effectively and efficiently processed by the preprocessing model, providing high-quality inputs to the prediction model. 3. The preprocessing step leads to a substantial increase in the accuracy of the prediction model. The prediction model based on NN has a much higher accuracy than the widely-used classical linear regression model in our 4D trajectory prediction task. A comparison study is conducted to demonstrate the effectiveness of the hybrid machine learning models. Comparisons are made on K-means++ and DBSCAN algorithms in the preprocessing part, Multiple Linear Regression (MLR) and NN model in the prediction part, as well as the combination of two parts.

Literature review

4D trajectory prediction has attracted numerous attention from worldwide researchers in the past decades. Considerable research efforts have been taken to enrich the 4D trajectory prediction approaches. From the view of methodology, the trajectory prediction models can be mainly classified into aircraft performance models and machine learning models [START_REF] Hrastovec | Prediction of aircraft performances based on data collected by air traffic control centers[END_REF].

Aircraft performance models belong to physics-based approaches. The model structure is based on kinematic assumptions. The model parameters are determined based on a model of the aircraft performance, the planned flight routes, the predicted atmospheric condition, and the expected command and control strategies given by pilots or Flight Management System (FMS) (known as Aircraft Intent). The most precise aircraft performance model is Base of Aircraft Data (BADA), which provides increased levels of precision in aircraft performance parameters for modelling and simulation [START_REF] Nuic | Bada: An advanced aircraft performance model for present and future atm systems[END_REF]. A variety of researches based on BADA and Aircraft Intent has been conducted. For example, [START_REF] Xi | Simulation study of algorithms for aircraft trajectory prediction based on ads-b technology[END_REF] presented a classified ADS-B-based trajectory prediction algorithm. Based on the state estimation by Kalman filter and intent information captured by a pretreatment and probability method, the aircraft trajectory can be predicted with computation efficiency and fewer errors. [START_REF] Porretta | Performance evaluation of a novel 4d trajectory prediction model for civil aircraft[END_REF] presented a novel aircraft performance model in consideration of the effect of wind, for aircraft lateral guidance and a new procedure for speed estimation. The model input includes navigation data and aircraft intent information, based on BADA set. Simulation results show that the model is suitable for reliable trajectory prediction. [START_REF] Kaneshige | Implementation of a trajectory prediction function for trajectory based operations[END_REF] described the implementation and evaluation of a motion-based trajectory prediction function, which can increase the resiliency and robustness of TBO. Based on the performance index, such as the fuel consumption, flight time, the algorithm computes the difference between with trajectory prediction and without trajectory prediction. Clearly, these aircraft performance models have made significant contributions to trajectory prediction, however, most of these models made ideal assumptions, rarely considering the actual constraints and human behaviour factors. The prediction in this method is conducted one by one, which cannot reflect the interaction between trajectories. In addition, once the data resources are limited or poorly supported, the prediction models will be ineffective for short-term trajectory prediction, not to mention large-scale trajectory prediction problems.

As a branch of Artificial Intelligence (AI), machine learning has been developed over 30 years. One of the major function of AI is learning from experience and make predictions. The trend of recent years shows that machine learning is commonly used in trajectory prediction domain. Compared with those aircraft performance models, machine learning models were built with weak assumptions or even without assumptions. In most cases, it shows better prediction performance. For example, Le [START_REF] Fablec | Using neural networks to predict aircraft trajectories[END_REF] used Neural Networks (NN) to predict an aircraft trajectory in the vertical profile. Trained by real historical trajectory dataset, 2 different model structures were presented. In the first method, the inputs include current altitude, the remaining altitude to reach Request Flight Level (RFL) and n past vertical speeds. The output is the next speed. The inputs of the second structure consist of starting altitude and remaining altitude to reach the RFL. n first initial speeds is the output. Simulation result showed that NN model gives better results than classical prediction functions based on the model of aircraft. [START_REF] Leege | A machine learning approach to trajectory prediction[END_REF] introduced Generalized Linear Model (GLM) for trajectory prediction at prediction horizons of 15NM to 45NM on fixed arrival route. The inputs of this model are aircraft type, ground speed at Initial Approach Fix (IAF), altitude over IAF, surface wind and altitude wind. All inputs come from surveillance data and meteorological data. [START_REF] Alligier | Machine learning applied to airspeed prediction during climb[END_REF] also applied GLM to predict the airspeed during climb. The simulation is conducted on actual Mode-C radar data of 9 different aircraft types. The prediction error is quite small, compared with BADA reference speed profiles.

In view of improving the accuracy in prediction tasks, [START_REF] Trivedi | The utility of clustering in prediction tasks[END_REF] carried out a study on the feasibility of utilizing clustering as a preprocessing approach. Their research demonstrates that the improvement on prediction accuracy is significant on large-scale clusterable datasets by combining the clustering with even some simple machine learning predictors. Under routine traffic situation in the TMA, the aircraft follows the standard arrival/departure procedure and regular ATC instructions, which makes trajectories clusterable. Thus, the application of machine learning together with clustering for 4D trajectory prediction is a valuable and interesting research topic. Numerous clustering models have been proposed for airspace monitoring, traffic patterns identification, automated traffic analysis, trajectory efficiency analysis [START_REF] Eckstein | Automated flight track taxonomy for measuring benefits from performance based navigation[END_REF][START_REF] Gariel | Trajectory clustering and an application to airspace monitoring[END_REF][START_REF] Chen | Clustering of trajectories based on hausdorff distance[END_REF][START_REF] Tang | A flight profile clustering method combining twed with k-means algorithm for 4d trajectory prediction[END_REF][START_REF] Ayhan | Diclerge: Divide-cluster-merge framework for clustering aircraft trajectories[END_REF][START_REF] Nicol | Unsupervised aircraft trajectories clustering: a minimum entropy approach[END_REF].

Based on these studies, models combining clustering with simple predictors have been investigated in the latest years. [START_REF] Tastambekov | Aircraft trajectory forecasting using local functional regression in sobolev space[END_REF] considered the short to mid-term aircraft prediction problem, namely, the prediction with a time horizon from 10 min to 30 min. The model firstly searches similar trajectories in terms of shape and time, then uses wavelet decomposition to solve the linear regression model of the relationship between time and trajectory projection onto one of the three axis X, Y and Z. This method gets good prediction results with high robustness and efficiency. [START_REF] Hong | Trajectory prediction for vectored area navigation arrivals[END_REF] introduced a new framework for predicting Estimated Time of Arrival (ETA) of aircraft from a specified entry fix. The model contains two steps: trajectory pattern identification by Dynamic Time Warping (DTW) and MLR models training for each pattern. The prediction of arrival time can be achieved by applying different MLR model for each trajectory pattern of target aircraft. However, most of the existed clustering-based prediction model are unable to conduct a comprehensive preprocessing step. Some models neglect the prediction steps, directly considering clustering results as prediction results. A majority of trajectory pattern identification approaches are not robust, requiring carefully selected flight data that follow the same departure/arrival procedure or pass the specified entry fix rather than in a specified airspace. If there exists some noise and overflights in the data, the predicted results will be far less effective. In addition, in the existed prediction model, the regression models that have been used are relatively simple and have a shallow structure. What's more, for different traffic patterns, the prediction processes have to be executed successively, without parallel computing ability.

Methodology

In order to improve the prediction accuracy and efficiency, we propose a novel hybrid trajectory prediction approach, see Fig. 1. There are two parts in this model: clustering-based preprocessing model and MCNN-based prediction model. The preprocessing step is mainly composed of re-sampling, data augmentation, PCA and DBSCAN algorithm. In this part, the model firstly cleans and reshapes the data, then identifies the 4D trajectories into different clusters and removes noise in an efficient way. Each cluster symbolizes that the corresponding trajectories have the similar pattern. Noise contains trajectories with holding patterns, trajectories with large vectoring, trajectories in special cases and overflight trajectories. Finally, the trajectory data quality will be highly increased. In the prediction step, we apply the MCNN model for 4D trajectory prediction. For each partition of trajectories, there is a predictor, in which there is an individual NN-based learning cell. Each individual learning cell will be trained with the associated cluster of trajectories. Consequently, each cluster of trajectories will have its corresponding prediction model. For the new input trajectory, we will classify them into different corresponding clusters by decision trees according to the entry point, then with already trained multi-cells predicting model, trajectory prediction of the input data will be generated.

Data preparation

We choose Automatic Dependent Surveillance-Broadcast (ADS-B) as the data source of the 4D trajectory prediction model. Each record of ADS-B data contains:

1. Type of operation (departure/arrival), 2. Coordinated Universal Time (UTC) time stamp t, 3. Aircraft number, 4. Position (X, Y, Z), 5. Heading Ψ, 6. Horizontal ground velocity V h , 7. Vertical ground velocity V v , etc.

Each record with the same aircraft number belongs to an aircraft i, and the collection of all records for that aircraft forms the trajectory T i , i = 1, ..., n, where n is the total number of trajectories in the dataset.

Clustering-based preprocessing

Data Cleaning and Formatting

To start with, the dataset is filtered on recording time scale, recording length and so on, making sure the data validity. Due to the instability of ADS-B receiver, the collected ADS-B data may not be continuous. Some trajectories may have missing parts. To solve this problem, it is necessary to place a low pass filter with the following function:

145 x1,i = x 1,i (1)
xl,i = αx l,i + (1 -α) xl-1,i , l = 2, ..., m i -1. ( 2
)
where the aircraft position in 3D coordinates and the heading of the l-th point of i-th trajectory are substituted into x l,i . α is a smoothing factor in [0, 1]. In this study, α is set to 0.5 to provide better results without too much delay. m i is the number of points in i-th trajectory.

In order to make dataset suitable for clustering, each trajectory should be represented as a vector. All the trajectory vectors are resampled into the same length, then their distance (which will be discussed in 150 the following section) can be computed. The resampling method for i-th trajectory is given as follows:

T i = T l,i l = floor k • m i S , k = 1, ..., S (3)
where the function floor(x) takes the largest integer that does not exceed x. l is a set containing all the index of sampled points for each trajectory, S is the number of sampled points for each trajectory. Experience shows that S = 50 would have been enough for clustering purpose. Trajectories with less than 50 points were eliminated due to statistical insufficiency.

Dimensionality augmentation

This step aims to augment the dimensionality of dataset for better clustering. The existing dimensions may not be sufficient and will result in lack of information, which can't completely reflect the differences between each trajectory. Thus, augmentation of dimensions will help improve the clustering performance. Therefore, the following dimensions will be added into the dataset:

R l,i = (X l,i -X ref ) 2 + (Y l,i -Y ref ) 2 + (Z l,i -Z ref ) 2 (4)
2. Distance from the corner point D. The corner point (X cor , Y cor , Z cor ) can be assigned as a corner of the considered airspace. The corner point will help solve the identifying problem when two trajectories are symmetric. The D l,i is calculated by the function below:

D l,i = (X l,i -X cor ) 2 + (Y l,i -Y cor ) 2 + (Z l,i -Z cor ) 2 (5)
The reference point and corner point play the role as multilateration. 3. Angular position from the reference point Θ. It shows the variation (turning status) of trajectories with respect to the reference point. Θ is defined as:

Θ l,i = arctan Y l,i -Y ref X l,i -X ref (6)
To sum up, the resampled dataset includes original features: UTC time stamp t, position (X, Y , Z), heading Ψ and additional features: distance from the reference point R, distance from the corner point D, angular position from the reference point Θ. To avoid the discontinuity at ±π and for feature scaling, the sine and cosine values of Θ and Ψ is utilised instead of the original value. Next, to make every feature on the same scale, each feature is normalised in [0, 1]. The general formula is given as:

x * = x -min(x) max(x) -min(x) (7)
where x is the original feature and x * is the normalised feature. Replacing x with our features, finally, all features of each trajectory is concatenated and organised as follows:

T i = [t * i P * i R * i D * i cos(Θ) * sin(Θ) * cos(Ψ) * sin(Ψ) * ] (8) T =    T 1 . . . T n    (9)
where

P * i = [X * i Y * i Z * i ].
Then, each trajectory is reshaped into a vector with 500 components. Matrix T is n × 500 and contains all trajectories with features, where n is the total number of trajectories.

Principal Component Analysis

As shown in Eq. ( 8), trajectories are related to various of factors. Nevertheless, among these factors, some are more related, while the other is less related. Redundant elements will decrease computational efficiency, even lead to larger errors. To solve this problem, PCA is introduced. PCA is a powerful tool used to reduce the dimension of dataset without losing too much information. The main idea of PCA is to derive an orthogonal linear transformation to project each of the vector variables into principal components for the maximum amount of variance that can be presented in lower dimensions [START_REF] Runkler | Data Analytics[END_REF].

PCA performs a linear transformation on the n × m (in this case m = 500) matrix T :

Y = E • T ( 10
)
where E is a rotation matrix, Y is the new principal component matrix. The variance of Y is:

var(Y ) = E T • C ( 11
)
where C is the covariance matrix of T , which can be written as:

C = 1 n -1 • T • T T (12)
The eigenvalues of C can be calculated as {λ i |i = 1, ..., m}, which correspond to the variances in Y as

{v i |i = 1, ..., m}, with λ 1 > λ 2 > ... > λ n .
To map a dataset X ⊂ R m to a dataset Y ⊂ R q with q ∈ {1, ..., m -1}, a rotation matrix E = (v 1 , ..., v q ) can be used. The dimension can be reduced by choosing the number of q. In order to keep most of the information while reducing dimensions effectively, it is required that the projection should better covers 90% of the variances, i.e., the cumulative percentage or variance explained G(q) is greater than 90%:

G(q) = q i=1 λ i m i=1 λ i ≥ 90% (13)

Clustering

Clustering algorithm belongs to unsupervised learning, for the purpose of partitioning unlabelled data according to similarity. The widely used clustering algorithm for air traffic pattern identification is K-means [START_REF] Eckstein | Automated flight track taxonomy for measuring benefits from performance based navigation[END_REF][START_REF] Tang | A flight profile clustering method combining twed with k-means algorithm for 4d trajectory prediction[END_REF] and DBSCAN [START_REF] Gariel | Trajectory clustering and an application to airspace monitoring[END_REF][START_REF] Chen | Clustering of trajectories based on hausdorff distance[END_REF], along with their improved or combination version.

The K-means algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] is an iterative clustering algorithm that partitions N observations into K clusters, where the number of clusters K is a hyperparameter. The K-means algorithm aims to minimise the Within-Cluster Sum of Squares (WCSS) cost function.

L(S; µ) = D K = K i=1 xj ∈Si x j -µ i 2 (14)
where S is the set of clusters partitioned by this algorithm, µ i is the mean of points in S i and K ≤ N . Initialisation of centroid seeds demands a-priori knowledge. The common idea is to use a farthest-first traversal approach or a selection at random from the dataset. Nevertheless, this is not robust to outliers. Instead, [START_REF] Arthur | k-means++: The advantages of careful seeding[END_REF] proposed K-means++ algorithm using a heuristic to initialise the centres. It improves both the computation speed and accuracy of K-means.

Density-Based Spatial Clustering of Applications with Noise (DBSCAN) is a commonly used densitybased clustering algorithm [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF]. The core concept of DBSCAN is to evaluate the density depending on the number of points within the -neighbourhood. DBSCAN classifies the observation points into three types: core point, density-reachable point and noise point. The algorithm expands to densityreachable areas from a selected core point, then obtaining a maximum area including the core point and density-reachable points. Being robust to the quality of datasets, DBSCAN can divide the dataset into several clusters and noise, where the a-priori selection of the number of clusters is not required. Besides, DBSCAN is able to find arbitrarily shaped clusters. The advantages of DBSCAN make it fit well with trajectory clustering scenarios. There are two hyperparameters in DBSCAN algorithm: the neighbourhood radius and the minimum number of points required to form a dense region minP ts. The value of affects the size of clusters. The value of minP ts affects the noise identification and the significance of clusters. After clustering, each cluster of trajectories can be considered as a main traffic pattern managed by controllers, and it has been passed by a large number of flights. The noise represents trajectories in some special cases, such as holding patterns and large vectoring.

In our case, the first q principal components of the organised trajectory matrix, computed by equation ( 13), are clustered by machine learning algorithm, DBSCAN or K-means++.

MCNN-based prediction model

In this paper, supervised learning techniques are used for short-term trajectory prediction. Each input x is a vector that contains the 3D position, heading, horizontal speed and vertical speed. The target variable t is the Estimated Time of Arrival (ETA) on the runway. We consider the following fixed-design regression model:

Y n = f (x n ) + n , n = 1, ..., N (15)
where Y n are the random variable that follows a mean function f (•) with errors n , which are independent and identically distributed (i.i.d.) random variables, such that E( n ) = 0 and var( n ) = σ 2 . The most classical machine learning model Multiple Linear Regression (MLR) aims to estimate the function f (•) under a linear assumption. Given multiple independent variables of N statistical units with D dimensions x n = (x n1 , ..., x nD ) T and corresponding dependent variable y n , n = 1, ..., N , the goal is to estimate the mean function f (•) by h w (•):

h(x; w) = D d=1 w d x nd + w 0 , n = 1, ..., N (16)
where {w n |n = 0, ..., N } are coefficients, which can be approximated by least squares regression, taking the estimator h w (•) that minimises the Mean Squared Error (MSE).

In our case, we use MCNN model to predict the ETA based on preprocessed real 4D trajectory data. For each cluster (cell) generated by the training set, an individual NN-based learning model will be built by training the dataset in the associated cluster. Consequently, each cluster will have an individual learning model. Then, the test set will be classified into different cell by classification algorithm. The output will be produced by the corresponding NN-based learning model. The advantage of using Neural Network (NN) in each cell is that they are able to learn the hidden and non-linear dependencies from the training data. The architecture of the proposed NN model for each cell is composed of an input layer, a hidden layer and an output layer, shown in Fig. 2. Given C cells, for each cell c, c ∈ {1, ..., C} and the input vector (x 1 , ..., x Nc ) T in test set T c , the output can be calculated as:

h c (x; w) = Ψ   Mc i=1 w i,c Φ   Nc j=1 w ij,c x j + w 0j,c   + w 0,c   (17)
where M c is the hidden layer node number, w ij,c is the weight between the j-th input node and the i-th hidden node, w i,c is the weight between the i-th hidden node and the output node, w 0j,c is the bias to the i-th hidden layer, w 0,c is the bias to the output layer, Φ is the sigmoid logistic function:

Φ(z) = 1 1 + e -z (18)
Ψ is the identity function:

Ψ(z) = z (19)
To reduce over-fitting, the error function we are trying to minimise is:

E(w) = (x,y)∈Tc (h c (x; w) -y) 2 + λ c Nc j=1 Mc i=1 w 2 ij,c (20)
where λ c is the regularisation coefficient for NN model in cell c. The back-propagation method was used to compute the gradient of the error function, and Levenberg-Marquardt algorithm [START_REF] Reynaldi | Backpropagation and levenberg-marquardt algorithm for training finite element neural network[END_REF] was used to minimise the error. The new input is classified by decision trees according to the entry point of the corresponding trajectory.

The advantages of the classification decision tree are that it does not acquire any a-priori knowledge and expensive computation costs. Once the classification tree model is built, the worst time complexity of classifying an unknown sample is Θ(H), where H is the maximum depth of the classification tree. In addition, in view of arrival flights in TMA, the longitude and latitude of entry points of trajectories in each cluster belong to a certain range in the airspace. This characteristic can be used to realise an effective classification on new input. In this work, Classification And Regression Trees (CART) algorithm is utilised for the decision tree, with the aircraft position X and Y as input.

Performance evaluation and model selection

To assess the performance of MCNN model, the MLR is proposed as the comparison learning model. DBSCAN and K-means++ are compared in preprocessing clustering part, so as to test different hybrid machine learning models.

In machine learning, hyperparameter stands for the parameters that are fixed before actual learning process begins, which cannot be directly adjusted by empirical risk. The bias-variance tradeoff is controlled by hyperparameters.

Firstly, in the preprocessing part, we will evaluate the performance of clustering algorithm DBSCAN, compared with K-means++. Different hyperparameters selection approaches are proposed for these clustering algorithms. The training data used is composed of the first q principal components of trajectory matrices T i in equation ( 8) that belongs to the training set of each outer cross validation fold. The number of clusters of K-means++ can be estimated by the silhouette method [START_REF] Rousseeuw | Silhouettes: a graphical aid to the interpretation and validation of cluster analysis[END_REF]. The silhouette value for each instance i is a measure of similarity between its cluster and of difference between other clusters. The silhouette value is defined as:

s i = a i -b i max(a i , b i ) (21)
where a i is the average distance from the i-th instance to the other instances of the same cluster, and b i is the minimum average distance from the i-th instance to instances in a different cluster. For each K in the hyperparameter grid, the value that maximises the average silhouette coefficient of all instances in the training set is chosen as the optimal cluster number of K-means++. On the contrary, the number of clusters of DBSCAN algorithm need not to be determined a-priori. However, DBSCAN algorithm has 2 hyperparameters and minP ts, in which can be evaluated by computing the K-Nearest Neighbour (K-NN) distance in the training set [START_REF] Sander | Density-based clustering in spatial databases: The algorithm gdbscan and its applications[END_REF]. The principal concept of this method is to calculate the distance of every point to its K nearest neighbour, and these distances are plotted in an ascending order.

The optimal value corresponds to the "knee" (threshold K-distance value), which is relatively simple to be recognised in a graphical representation. Since K-NN is a non-parametric method, the value of K is specified manually. With regards to minP ts, there is no general approach to estimate the value. In order to remove most abnormal trajectories and balance the size of dataset in the meantime, we search an optimal value of minP ts in the hyperparameter grid to make sure that the noise account for 5% ∼ 8% of the dataset. Secondly, in order to well select the hyperparameters of prediction model, and to achieve an unbiased performance of the prediction model, this paper uses nested cross validation method. It consists of the outer loop and the inner loop. A K 1 -fold cross validation splits the dataset S into K 1 subsets S i , i ∈ {1, ..., K 1 }.

For each outer iteration i, K 1 -1 folds S -i = S\S i act as training sets and one fold S i is test set. Then, there is another K 2 -fold cross validation, which will further split the training sets S -i into K 2 subsets S -i,j , j ∈ {1, ..., K 2 }. For each inner iteration j, K 2 -1 folds S -i \S -i,j play the part of training sets and the remaining fold S -i,j is validation set. The purpose of the inner loop is the selection of hyperparameters and the outer loop aims to assess the model performance. Taking K 1 = 5, K 2 = 5, the proportion of training sets, validation sets and test sets is 64%/16%/20%. The kernel idea of the whole process can be summarised on Fig. 3.

Here, we use a grid search hyperparameters tuning algorithm [START_REF] Alligier | Machine learning applied to airspeed prediction during climb[END_REF] on NN-based learning model. This algorithm aims at selecting the hyperparameter λ of an algorithm A λ by performing a 5-fold cross-validation on a set of examples T , which is depicted by algorithm 1.

Algorithm 1 Hyperparameters tuning

1: function TuneGrid(A λ , grid)[T ] 2: λ * ← arg min λ∈grid CV 5 (A λ , T ) 3: return A λ * [T ] 4: end function
The coefficient of determination R 2 and the statistical significance p-value are considered as the statistics for the error between the actual and estimated results of the MLR model. R 2 is defined as equation ( 22), providing a measure of how close the dependent variables are to the independent variables. The p-value tests the null hypothesis that the coefficient is equal to zero.

R 2 = 1 - N i=1 (y n -ŷn ) 2 N i=1 (y n -ȳ) 2 (22)
The performance of classification model is estimated by K-fold cross validation on the training set, setting K = 5 in this case.

Finally, we use Mean Absolute Error (MAE) and Root Mean Squared Error (RMSE) to assess the

MAE = 1 n n i=1 |ŷ i -y i | (23) RMSE = 1 n n i=1 (ŷ i -y i ) 2 , ( 24
)
where ŷi is the i-th predicted value and y i is the i-th observed value of ETA. Given that each outer iteration produces a MAE i , and a RMSE i , i = 1, ..., k 1 , the average MAE and RMSE can be computed as follows:

MAE = 1 k 1 k1 i=1 MAE i , (25)
RMSE = 1 k 1 k1 i=1 (RMSE i ) 2 , (26)

Experiment results and discussion

Dataset

The available dataset includes ADS-B records in July, 2017 over the TMA of Beijing Capital International Airport (BCIA), which is one of the busiest airport in the world, with three parallel runways: 18R/36L, 18L/36R and 01/19. Since the studied airspace is relatively small, the longitude, latitude of trajectory points in WGS84 can be projected into a projection plane. Note that, in this paper, arrival flights that correspond to QFU 36 are taken into consideration, where QFU means the magnetic orientation of runway-in-use. QFU 36 goes to North.

The dataset that we used in the experiments contains 8,677 arrival flights of QFU, being made up of 902,457 trajectory points. After data cleaning procedure, 103 flights were removed from the dataset, and the data validity rate is higher than 98%. Then, the dataset is formatted and augmented with dimensionality, depicted in section 3.2.

Clustering performance

First, in this study, the cumulative percentage of variance of trajectory matrix in equation 8 is calculated and presented in Fig. 4a. When the principal component reaches over 33, the variance explained will exceed 90%. Let q = 33, then the dimension of each trajectory was reduced to 33 from 500. To sum up, dimensionality augmentation enriches the features that principal components can choose. PCA reduces the dimension of the dataset, which makes the following clustering step more efficient and accurate in the projected principal component space. Then, the hyperparameters grids of aforementioned machine learning approaches are used, see Table 1. According to the density distribution of the dataset, K is set to 5 in K-NN. In all outer cross validation folds, K-NN distance plots of instances in the training set indicate that = 6.5. Fig. 4b represents the K-NN distance plot in one of the folds. minP ts = 50 satisfies that the noise accounts for 5% ∼ 8% of the dataset. Among the hyperparameter grid of K-means++, K = 4 minimises the average silhouette value of all instances in the training set for every outer cross validation. Fig. 4c is extracted from one of the cross validation fold as an example.

We find that DBSCAN algorithm clusters the trajectories into 5 partitions and noise. Taken one of cross validation folds into consideration, Fig. 5a represents 3D plot for the first principal components of the trajectory matrix. In this case, cluster points involve 93.00% and the noise points account for 7.00%. Black crosses indicate noise, and circles with different colour stand for trajectories in different clusters. Fig. 5b demonstrates the corresponding 2D trajectories. It can be seen that 5 main arrival procedures are clustered. Fig. 5c shows trajectories recognised as noise, which is mainly composed of holding patterns and trajectories with large vectoring, which will have an interference for the training stage. Therefore, the noise should be deleted from the dataset. With the same cross validation fold as Fig. 5, K-means++ clusters all the trajectories into 4 partitions. Unlike DBSCAN, K-means++ would not identify the noise. Fig. 6a represents the first 3 principal components of the trajectory matrix. Fig. 6b shows the corresponding trajectories, in which the cluster #4 combine cluster #3 and cluster #5 in Fig. 5b together. The clustering result of K-means++ is not good, considering the combining cluster contains several arrival procedures.

In addition, Fig. 7 the shortest transit time in TMA, which is from the west with less deviation before landing. Besides, in Tab. 3, compared with DBSCAN-based models, we find that the average number of nodes of classification tree for K-means++ is 11 times larger, the classification accuracy for K-means++ is about 3% lower, and the percentage variance of clusters in different outer cross validation folds of K-means++ is much less stable. Taking Fig. 5 as an example to explain the reason of such differences. In Fig. 5a and Fig. 5b, there are some clusters that intersect. In addition, because K-means++ cannot identify the noise, some abnormal or deviated trajectories were partitioned in clusters. Therefore, the classification trees have to devote more efforts for noise and clusters with intersection, which need more binary classifiers for fine classification that may cause overfitting. On the contrary, the shape of clusters generated by DBSCAN is quite good, and there is no obvious intersection between clusters. Even the decision trees for K-means++ are more complicated, the overall performance is getting worse.

Conclusion

In this paper, a novel 4D trajectory prediction approach that combines clustering model with prediction model is proposed, implemented and simulated for ETA prediction. The proposed model contains clusteringbased preprocessing step and MCNN-based machine learning prediction step.

First of all, four steps in the preprocessing part are discussed in detail, including data cleaning and Numerical experiments demonstrate that NN with DBSCAN preprocessing performs best in terms of MAE and RMSE, compared with other 5 methods, including NN with K-Means preprocessing, NN without preprocessing, MLR without preprocessing, MLR with K-Means preprocessing and MLR with DBSCAN preprocessing. Clearly, the hybrid machine learning prediction model MCNN with DBSCAN preprocessing can make an accurate, efficient and robust short term trajectory prediction in TMA.

Future work could be conducted on a larger time scale or special scale. Moreover, limited by the size of dataset that we obtain at this moment, deep learning models are not used in this paper. More advanced

Figure 1 :

 1 Figure 1: Proposed 4D trajectory prediction approach

  1. Distance from the reference point R, which indicates the convergence of aircraft toward the center of the TMA. Due to the runway configuration, we can define the reference point (X ref , Y ref , Z ref ) as the centre of the runway. For each trajectory point, R l,i is given as:

Figure 2 :

 2 Figure 2: Neural network architecture used in this paper

Figure 3 :

 3 Figure 3: Nested cross validation procedure

Figure 4 :

 4 Figure 4: Parameters selection of algorithms in the preprocessing step

  illustrates the decision trees for classification according to clusters generated by DBSCAN and K-means++ in one of the outer cross validation folds. It is clear that according to maximum depth and the number of nodes, the time complexity and the overfitting risk of Fig. 7a are much smaller than Fig. 7b.

Figure 5 :F

 5 Figure 5: DBSCAN clustering example in one of the outer cross validation folds

Figure 6 :Figure 8 :

 68 Figure 6: K-means++ clustering example in one of the outer cross validation folds

Table 1 :

 1 The grids of hyperparameters for machine learning algorithms

	Method	Hyperparameter grid
	K-means++	K = {2, 3, 4, 5, 6, 7, 8, 9, 10}
	DBSCAN	minP ts = {10, 20, 30, 40, 50, 60, 70, 80, 90, 100}
	MCNN	λc = {0.0001, 0.001, 0.01, 0.1, 0.5, 1, 2, 5}, c = 1, ..., C
		Mc = {2, 3, 4, 5, 6, 7, 8, 9, 10}, c = 1, ..., C

  Compared with the performance of ETA prediction on NN and MLR with K-means++-based preprocessing step on each outer cross validation fold, It is also indicated that NN model performs still better than MLR model for each cluster in terms of MAE/RMSE. Last but not least, in Tab. 4, regarding overall prediction errors MAE and RMSE, the NN model with DBSCAN preprocessing performs best, followed by NN model with K-mean++. MLR model without preprocessing performs worst on the prediction task.

Table 2 :

 2 The performance of ETA prediction on NN and MLR with DBSCAN-based preprocessing step on each outer cross validation fold

	Outer CV#	Classification tree	Cluster#	Percentage	DBSCAN P. + NN	DBSCAN P. + MLR
		Accuracy	Size			MAE[s]	RMSE[s]	MAE[s]	RMSE[s]	R 2	p-value
				1	5.15%	134.20	166.56	162.80	201.44	0.90	≈ 0
				2	12.19%	192.25	275.28	229.83	314.99	0.90	≈ 0
	1	99.92%	19	3	54.98%	128.60	184.52	154.94	215.48	0.88	≈ 0
				4	11.71%	84.31	142.19	90.32	144.31	0.94	≈ 0
				5	9.26%	161.13	223.92	200.39	255.23	0.89	≈ 0
				1	5.42%	164.42	218.15	175.22	236.19	0.90	≈ 0
				2	12.03%	184.90	286.96	220.73	312.54	0.91	≈ 0
	2	99.61%	15	3	55.74%	142.85	233.94	165.28	251.08	0.88	≈ 0
				4	11.21%	106.91	160.81	120.16	173.04	0.94	≈ 0
				5	9.48%	188.87	356.39	218.93	368.60	0.92	≈ 0
				1	4.93%	128.23	161.03	134.21	176.84	0.90	≈ 0
				2	10.98%	161.44	203.57	236.30	287.21	0.91	≈ 0
	3	99.84%	19	3	56.66%	144.49	211.19	162.70	228.78	0.88	≈ 0
				4	10.76%	81.07	108.45	113.08	144.20	0.94	≈ 0
				5	9.40%	146.32	356.39	195.29	249.01	0.89	≈ 0
				1	4.99%	122.03	155.63	153.30	202.44	0.90	≈ 0
				2	11.12%	160.28	218.39	198.74	268.66	0.91	≈ 0
	4	99.77%	21	3	57.59%	123.23	173.47	137.96	183.15	0.88	≈ 0
				4	10.50%	95.19	151.39	154.52	209.90	0.94	≈ 0
				5	9.86%	139.49	293.05	185.23	330.03	0.88	≈ 0
				1	5.02%	124.38	155.63	147.56	204.50	0.90	≈ 0
				2	11.72%	157.05	227.16	193.07	252.01	0.91	≈ 0
	5	99.92%	11	3	56.01%	123.80	182.80	143.88	200.72	0.87	≈ 0
				4	10.86%	117.85	333.58	165.29	355.82	0.94	≈ 0
				5	9.39%	144.61	183.42	171.50	220.41	0.88	≈ 0

Table 3 :

 3 The performance of ETA prediction on NN and MLR with K-means++-based preprocessing step on each outer cross validation fold

	Outer CV#	Classification tree	Cluster#	Percentage	K-means++ P. + NN	K-means++ P. + MLR
		Accuracy	Size			MAE[s]	RMSE[s]	MAE[s]	RMSE[s]	R 2	p-value
				1	56.58%	132.95	181.09	149.10	205.64	0.86	≈ 0
				2	13.00%	182.60	264.28	216.32	297.36	0.88	≈ 0
	1	96.57%	189	3	18.21%	119.65	189.74	124.61	173.45	0.88	≈ 0
				4	12.20%	167.83	248.05	201.59	279.07	0.88	≈ 0
				1	58.07%	141.21	235.99	158.17	246.24	0.86	≈ 0
				2	11.18%	184.53	345.09	217.16	348.63	0.92	≈ 0
	2	98.61%	121	3	17.85%	128.03	186.80	133.91	191.13	0.90	≈ 0
				4	12.90%	187.92	286.01	217.12	306.47	0.89	≈ 0
				1	57.99%	141.89	209.23	166.78	233.82	0.86	≈ 0
				2	11.98%	164.78	210.26	250.49	301.44	0.88	≈ 0
	3	95.85%	225	3	17.13%	102.30	133.47	127.53	160.14	0.87	≈ 0
				4	12.90%	157.33	209.05	201.26	261.95	0.88	≈ 0
				1	58.36%	118.65	162.36	141.73	187.14	0.85	≈ 0
				2	12.09%	161.72	214.28	198.10	268.12	0.89	≈ 0
	4	96.79%	203	3	12.80%	115.25	146.73	127.20	164.58	0.86	≈ 0
				4	16.75%	131.08	230.94	137.81	242.00	0.87	≈ 0
				1	45.36%	125.98	191.33	138.49	201.65	0.86	≈ 0
				2	17.28%	135.32	232.01	140.62	226.50	0.88	≈ 0
	5	96.06%	195	3	24.68%	150.65	243.84	157.07	257.20	0.88	≈ 0
				4	12.68%	162.09	226.63	190.85	253.26	0.88	≈ 0

Table 4 :

 4 The average performance of ETA prediction on NN and MLR with/without preprocessing step over 5 folds

	Preprocessing	Classification tree	NN		MLR
		Accuracy	Avg number of nodes	MAE[s]	RMSE[s]	MAE[s]	RMSE[s]	R 2	p-value
	Without	-	-	145.58	218.26	222.88	302.59	0.77	≈ 0
	K-means++	96.78%	187	139.09	212.20	160.01	232.76	0.81	≈ 0
	DBSCAN	99.81%	17	135.62	210.21	162.85	236.04	0.81	≈ 0

model prediction performance at different time horizons. In each sub-figure, the X-axis presents the transit time t(s) in TMA. While t(s) = 0, it means that aircraft is entering the TMA, and t(s) = 800 means that aircraft is flying in TMA for 800 s. Note that the shortest transit time in TMA is around 800 s, which mainly corresponds to the flights entering from west. When t(s) > 800, some of the flights may have landed on the runway. Fig. 8 shows the prediction errors for trajectories from all the entry points of TMA. The Y-axis depicts the absolute prediction error between predicted ETA and observed Actual Time of Arrival (ATA). For each 50 s time step, a box-plot shows the 5%, 25%, 50%, 75% and 95% quantiles.

The remaining 10% error mostly come from abnormal trajectories (highly deviated, go around, with holding procedure, in emergency, etc.), which are difficult to be predicted. So this part of error is not displayed. Two models are presented in each sub-figure. For clearly visualisation, two models are overlaid. We find that with different prediction methods, the results show the same trend, that is: when the time to destination is fewer, the absolute prediction error is smaller. We further observe the distribution of ETA prediction errors. Note that the whiskers contains 90% of the prediction errors, and the box contains 50% of the prediction errors. At the entry time of For different traffic patterns (clusters), the prediction error from different entry point of TMA (t(s) = 0) is different. The statistics of prediction results on NN and MLR with DBSCAN-based preprocessing model in each clusters and each outer cross validation fold were summarised in Tab. 2 and Tab.4. For MLR model, the correlation coefficient R 2 between the ATA and estimated ETA for each clusters is computed, representing the percentage of the variance explained in the ATA. The computed p-value for the null hypothesis was almost zero, which was small enough to reject the null hypothesis at a significance level of 0.01.

In Tab. 2, new test data can be well clustered into 5 partitions by application of the classification tree and the DBSCAN algorithm. The accuracy of the classification tree is more than 99.5%, and in view of MAE and RMSE, the proposed NN model performs better than the MLR model for each cluster. For cluster #4, the associated MAE is generally smaller than the other partitions, on account of the shape and the flying distance of the arrival procedure of trajectories. We also find that cluster #4 usually has deep learning model can be applied to improve the prediction performance when the dataset is enough for training the model.

Acknowledgement

The authors would like to thank Serge Roux who provided us with the data used in the computational experiments. Special thanks also go to colleagues in the Optim group of ENAC for their generous suggestions.