Unstable Approach Detection and Analysis Based on Energy Management and a Deep Neural Network

Abstract

The study of managing risk in aviation is the key to improving flight safety. Compared to the other flight operation phases, the approach and landing phases are more critical and dangerous. This study aims to detect and analyze unstable approaches in Taiwan through historical flight data. In addition to weather factors such as low visibility and crosswinds, human factors also account for a large part of the risk. From the accidents studied in the stochastic report of the Flight Safety Foundation, nearly 70% of the accidents occurred during the approach and landing phases, which were caused by improper control of aircraft energy. Since the information of the flight data recorder (FDR) is regarded as the airline’s confidential information, this study calculates the aircraft’s energy-related metrics and investigates the influence of non-weather-related factors on unstable approaches through a publicly available source, automatic dependent surveillance-broadcast (ADS-B) flight data. To evaluate the influence of weather- and non-weather-related factors, the outliers of each group classified by weather labels are detected and eliminated from the analysis by applying hierarchical density-based spatial clustering of applications with noise (HDBSCAN), which is utilized for detecting abnormal flights that are spatial anomalies. The deep learning method was adopted to detect and predict unstable arrival flights landing at Taipei Songshan Airport. The accuracy of the prediction for the normalized total energy and trajectory deviation of all flights is 85.15% and 82.11%, respectively. The results show that in different kinds of weather conditions, or not considering the weather, the models have similar good performance. The input features were analyzed after the model was obtained, and the flights detected as abnormal are discussed.

1. Introduction

With the increase in global air traffic, aviation safety issues have come into focus in various countries. The Global Aviation Safety Plan 2020–2022 edition proposed by the International Civil Aviation Organization (ICAO) [1] has recognized five high-risk categories (HRCs) as global priorities, namely controlled flight into terrain (CFIT), loss of control in-flight (LOC-I), mid-air collision (MAC), runway excursion (RE), and runway incursion (RI). Among these five high-risk categories, CFIT, LOC-I, and RE are the most common types of approach and landing accidents. According to statistics from the ICAO safety report 2021 edition [2], accidents caused by RE made up 18.8 percent of all accidents and 29 percent of accidents with planes mostly damaged or destroyed in 2020. These accidents with 24 fatalities also represented 50 percent of all fatal accidents. There is one accident associated with LOC-I which represented 2.1 percent of the total accidents with zero fatalities. According to Taiwan Aviation Occurrence Statistics (2010–2019) [3], which is provided by the Civil Aeronautics Administration (CAA), the flight phase in which most events occur is the approach and landing phase. Since pilots need to make lots of decisions and apply many actions during the approach and landing phase, the most common factors causing the events were human-related, which stands at 46.7 percent (35 out of 75 accidents). The mandatory report provided by the CAA is a report that records the events that have occurred at Taiwan’s airports or were caused by Taiwanese airlines around the world. It tends to figure out the causation of the factors and improve flight safety. However, the mandatory report does not provide enough information to investigate the factors, which makes the related causes of the events that occur at Taiwan’s airports ambiguous. Therefore, this study intends to find out the relationship between the events and the factors and the interaction between the weather- and non-weather-related factors, which are the flight trajectories and energy-related metrics of the flights.

According to Flight Safety Foundation ALAR Briefing Note 4.2 [4], either inadequate energy or excess energy may cause an approach and landing incident or an accident involving loss of control, abnormal landing, tail strike, or runway excursion. There are statistical data showing that approximately 70 percent of the approach and landing accidents that had few survivors were caused by inadequate energy management. According to the International Air Transport Association (IATA) safety report 2021 edition [5], the contributing factor that has the highest proportion in the flight crew errors category is manual handling/flight controls, which exerted an influence on 38 percent of aircraft accidents that happened in 2021. Research by the U.S. National Transportation Safety Board (NTSB) [6] stated that a high airspeed at the final approach fix (FAF) may prevent the capture of the specified gliding path by the autopilot and may impede the aircraft from stabilizing at the expected stabilization height.

To improve flight safety, information technology has been applied to process and analyze the recorded or received flight data. Flight operational quality assurance (FOQA), known as flight data monitoring (FDM), is the most used method which works by capturing and analyzing flight data to find solutions to improve flight safety and increase overall operational efficiency. As a result of an ICAO Annex 6 mandate, all airlines are required under regional legislation to implement FDM programs. However, FOQA analysis is based on accurate flight data obtained from the flight data recorder (QAR), but the information of the FDR is regarded as the confidential property of airlines. Fortunately, the publicly accessible ADS-B data are an alternative source of flight information and have been widely used in many applications, such as anomaly detection, air traffic flow identification, and trajectory analysis. Aircraft obtain accurate flight information including position and velocity data via satellite navigation or other sensors and periodically broadcast it via ADS-B. Even if ADS-B contains less information than FDR data, its accessibility and continuity in data offering provide it with more opportunities to carry out more research.

In the past, most studies on flight safety have focused on the detection of outlier flights or taken the characteristics of accidents to identify deficiencies. Along with improvements in data collection and the industry moving to carry out a more proactive approach to improve flight safety, machine learning algorithms, or deep learning, have come into sight in data analysis and provided more prospective information such as the forecast of flight trajectory or the prediction of flight parameters. Many works have analyzed the safety of the approach and landing phases with unsupervised learning methods such as clustering frameworks, density-based spatial clustering of applications with noise (DBSCAN), and SVM [7,8,9,10,11]. Corrado [10] analyzed anomalies through spatial and energy aspects by applying HDBSCAN and DBSCAN and discussed the relationship between spatial anomalies and energy anomalies. They had good results in detecting anomalous flights and abnormal flight parameters. However, there is no accredited standard for the threshold of identifying normal and abnormal in unsupervised learning, and it is therefore decided by the author, who may exert subjective judgment.

The previous work of this study was to develop a supervised deep learning method, analyzing the risk of the final approach phase caused by weather and determining a risk index for the landing flight [12]. The flight path angle during the approach and landing phases is an indicator of the risk assessment, which is an output of the deep neural network model, while the inputs are the weather conditions and the trajectory deviation. The author displayed a nice predicative accuracy of the flight path angle. Nevertheless, since weather is just one of the factors influencing approach safety and there are more important factors that interfere that have not been considered, the safety assessment of the approaching phase can be discussed one step further. The pilot’s operations or their reaction when facing poor weather is the main reason that puts flight safety at risk. As a result, people have gradually attached importance to research on the energy management of aircraft, which will be discussed in the later paragraphs.

At present, most of the studies related to factor analysis of unstable approaches use exceedance analysis to identify the events and factors involved. However, exceedance analysis cannot consider human–automation interaction problems, and it is difficult to explain the relationship between the events and the factors [13]. Energy-based metrics such as potential energy and kinetic energy, which quantify the aircraft’s energy state and can be regarded as an objective factor for evaluating various safety-critical conditions, will be suitable tools to identify unstable approach events and to explain the connection between incidents and causes. Consequently, the appropriate energy management of an aircraft is a critical subject for enhancing aviation safety. The flight crew’s incapability to evaluate or supervise the aircraft’s energy state during the approach phase is often cited as a reason for unstable approaches. A significant number of fatal accidents have been associated with poor management of vertical flight paths (potential energy) and/or airspeed (kinetic energy) [14] Airbus [15] has pointed out that energy management is the key issue to be discussed and would play an important role in improving flight safety.

Sembiring [16] presented a generic approach using QAR data for detecting unstable approaches through energy management which can be used for benchmarking purposes. Since every airline has its own set of rules for detecting unstable operations, the rule-based method cannot be a united standard. The author defined an unstable approach detection framework starting from the initial approach point, which is about 10 nautical miles away from the touchdown point to the touchdown point, with the concept of the minimum stabilization height under instrument meteorological conditions (IMC). The standard energy state is a key to energy management, so the author constructed a structure with the ideal landing energy curve and upper and lower energy boundaries. The formula of total energy, which equals the addition of potential energy and kinetic energy, was used to define the energy curves and identify unstable approaches through energy boundary construction. A quantitative feature of the curve would be a measurement of detecting unstable approaches, whereby the excess area outside of the boundary was chosen by the author as the classification rule. When the excess area of a flight energy curve exceeds a threshold value, it would be detected as an unstable approach. This unstable detection method based on energy management required fewer than ten recorded parameters, which is valuable when there is little information, and led to a generic formulation for detecting genuine unstable approaches.

Puranik [7] and Kumar [17] used energy-based metrics with the unsupervised machine learning method for identifying anomalous flight data records and trajectories in the approach and landing phases, with the use of FDR data and ADS-B data, respectively. Since energy-based metrics make the aircraft’s energy state measurable, they can be viewed as an objective indicator to evaluate diverse safety-critical situations. Using energy metrics, specific features that deviate outside the normal boundary may be individually selected for further inspection. Four different subsets of the input feature vector were compared with each other in Puranik’s study [7], with them detecting anomalies in 1000 flights. The results showed that even with a limited amount of information available and using different energy-metric subsets for evaluating energy metrics, the methodology can identify 80% of anomalies. In Kumar’s study [17], the author identified air traffic flows with the HDBSCAN method before energy anomaly detection. Through the isolation forest algorithm, anomalies were detected, and an anomaly score was given to better understand the severity and the relative level of anomaly compared to other trajectories.

From the above introduction, there is no method to identify the factors influencing events; therefore, this study intends to find out the factors that are non-weather related, since human-related factors need more precise data which are hard to obtain. The objectives of this dissertation are as follows:

• Identify whether the unstable approach factors are weather-related factors or non-weather-related factors.
• To detect the unstable approach through energy metrics and trajectory deviation.

The remainder of this study is organized as follows. Section 2 presents the research process and an overview of the data used, and the methodologies applied in this paper are also presented. In Section 3, the details of data preprocessing and the preparation works before the deep learning model training process are shown. The introduction of the deep neural network, the model architecture, the model tuning process, the feature importance analysis, and a discussion of the analysis results are shown in Section 4. The conclusions of this study and suggestions for future work are addressed in Section 5.

2. Research Process and Data Analysis

This paper aims to identify unstable approaches by calculating the flight’s energy metrics and finding out the factors that cause unstable approaches. By using the energy concept for identifying unstable approaches, the airport condition will not be considered, and thus, the method can be used at all airports for ATR 72-600 landing flights in Taiwan. Section 2.1 introduces the research process of the whole study. After that, Section 2.2 introduces the data used in this paper.

2.1. Research Process

The research flowchart of this study is shown in Figure 1. The ADS-B data from the website FlightRadar24 are the main data used in this research. After the data were extracted and cleaned, the corresponding weather condition of each flight could be found in the meteorological terminal aviation routine weather report (METAR). The weather condition of each flight was labeled by applying the ATMAP classification algorithm [12], and the flights landing on RWY 10 at Taipei Songshan Airport were classified according to the weather label of each weather condition. The outliers of each group classified by weather labels were detected and eliminated from the analysis by applying HDBSCAN afterward. Then, the flight parameters of each flight, which are the input features of the model, were calculated by using ADS-B data (altitude, speed, vertical speed, etc.). The trajectory deviation of each flight was calculated by the position data from ADS-B, which is one of the model’s outputs. Another output of the model is the normalized total energy of the flight, which is computed by the aircraft’s normal and ideal landing speed, V_ref. Afterward, the deep neural network can be trained, the feature importance can be calculated, and the model will be validated by the real events as the last step.

2.2. Data Collection and Analysis

2.2.1. Data Sources and Description

The ADS-B data of the ATR 72-600 aircraft landing at Taipei Songshan Airport (IATA: TSA; ICAO: RCSS), which originated from the website FlightRadar24, were used in this study for model training, and the time interval is a year (1 July 2019~30 June 2020). The timestamp, UTC, callsign, position (latitude and longitude), altitude, speed, and direction are the attributes involved in ADS-B data, as shown in

Table 1. ADS-B attributes description.

Attributes	Description
Timestamp	A sequence of characters or encoded information, usually giving the date and time of day
UTC	Coordinated universal time
Callsign	Consists of the telephony designator of the aircraft operating agency, followed by the flight identification
Position	Latitude and longitude of the aircraft
Altitude	Calibrated altitude
Speed	Ground speed
Direction	Aircraft’s heading, expressed in degrees from true north

. Latitude, longitude, and altitude will be used for trajectory deviation calculation. Information such as speed and altitude will be used for defining the aircraft’s energy state, where speed is essential for kinetic energy and altitude is necessary for potential energy.

2.2.2. Weather Data Labeled by the ATM Airport Performance (ATMAP) Weather Algorithm

The landing flights will be classified according to similar weather conditions beforehand to exclude the influence of weather when finding the factors affecting the approach. The weather data from the METAR, which are standardized by the ICAO to be in common use around the world, are labeled by the ATMAP weather algorithm [12]. The ATMAP weather algorithm is a unified evaluation of weather conditions at European airports [18]. Weather phenomena, weather class, severity code, and coefficient are the basic elements of ATMAP weather algorithms, as shown in Figure 2. A severity code with a coefficient is assigned to each weather class under observation. Nevertheless, the ATMAP weather algorithm was designed for the European area, so the measurement does not perfectly reflect the weather situation in Asia, so a part of the classification criteria must be modified to better fit the situation for Taiwan’s airports [12].

Table 2. Weather phenomena and descriptions.

Weather Phenomena	Descriptions
Visibility	A measure of the opacity of the atmosphere at Taipei Songshan Airport
Ceiling Height	The height above the Earth’s surface of the lowest layer
Wind Speed	Wind observed in local routine reports used for arriving or departing aircraft
Crosswind	The wind direction that has a perpendicular component to the direction of travel
Precipitation	Any of the forms of water particles, whether liquid or solid
Temperature	Temperature
Dangerous Phenomena	A harsh weather condition
CB/TCU without precipitation	Cumulonimbus clouds and towering cumulus without rain
CB/TCU with precipitations	Cumulonimbus clouds and towering cumulus with rain
Included Angle	The angle between the wind direction and the landing runway

. shows the weather conditions that were obtained from the METAR report of Taipei Songshan Airport and considered in this study.

2.2.3. Airport Information and Navigation Procedures

Taipei Songshan Airport has scheduled flights to domestic destinations in Taiwan and a small number of international flights. There is one runway at Songshan Airport, with a length of 2605 m, and the directions of the runway are 100 (designator: 10) and 280 (designator: 28), with the runway threshold elevations being 13 ft and 18 ft, respectively. The number of aircraft movements was about sixty thousand in 2019 but decreased rapidly afterward, especially international flights, down to forty thousand in 2020 due to the influence of the COVID-19 pandemic. The distribution of the aircraft movement of the TSA includes domestic flights, which account for about ninety percent, and international and cross-strait flights, which account for ten percent of the aircraft movements.

2.3. Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN)

Clustering based on density is a common way to group data (e.g., DBSCAN, OPTICS, and DENCLUE). For estimating the density of a dataset, a parameter needs to be decided: the core distance, the distance of a point to its nearest neighbor. In denser regions, the core distance will be smaller than in the sparser region, and it will be big for the core distance. The core distance is what makes these methods “density-based”.

However, these existing methods have several limitations. Firstly, some methods can only provide non-hierarchical labeling of the data objects, based on a global density threshold. A single-density threshold for clustering cannot suitably depict datasets with clusters of various densities and/or nested clusters. Secondly, many hierarchical methods only imply how to extract a flat partition using a density threshold. The most important clusters may not arise if these clusters are in different density flats [19]. The dataset that is going to be clustered is usually not separated evenly, its shape may be arbitrary, with different sizes or densities, and there might be outliers. A robustness algorithm with fewer assumptions is required for real-world data exploration.

If there is only one density threshold for clustering, the method will either over-group the clusters with different densities visually or will fail to include the entire cluster with identical density. HDBSCAN first builds a hierarchy to find out which densities of the clusters merge last and in what order, and then for each cluster, HDBSCAN can decide whether to divide a cluster into its subclusters. Through these procedures, it automatically generates the clusters. Besides clustering, HDBSCAN can detect the outliers in the dataset. Not only will the points at the outlying border be highlighted as outliers but so will the points between the clusters if the core parts of the clusters are far enough. The primary parameter that needs to be decided when using HDBSCAN is the minimum cluster size (min_cluster_size), which is the smallest quantity of data to form a group (cluster).

2.4. Energy Boundary Construction

Since the management of the aircraft’s total energy during the approach affects the stability, it is necessary to focus on the energy factors of the aircraft. In this study, the total energy composed of kinetic energy and potential energy is one of the indicators of identifying an unstable approach [16]. The kinetic energy of the aircraft is determined by the aircraft’s speed, and the potential energy is decided by the aircraft’s altitude or the aircraft’s gliding angle. The energy used in energy management is based on energy height, which is also termed specific energy in the study [10] and defined as real energy divided by gravity and has the same unit as height [16]. The specific total energy is defined as follows:

$$\mathrm{Specific}\mathrm{Total}\mathrm{Energy}=h_{AGL}+\frac{V^2}{2g}$$

(1)

where $h_{AGL}$ is the height above ground level, V is the aircraft’s speed, and g is the gravity constant.

2.4.1. Energy Boundary

The energy boundary is constructed by defining the upper and lower boundary of the total energy, and the ideal status of the aircraft’s total energy is used for boundary normalization [16]. The upper energy boundary of the aircraft is the maximum total energy that the aircraft can have during the approach phase. The definition of the upper boundary is different when the aircraft is below and above the IMC gate. The IMC gate is an important stabilization checkpoint under instrumental meteorological conditions, which is 1000 feet above ground level for aircraft with an ILS. When the aircraft is below the IMC gate, the speed of the aircraft, which is the kinetic part of the total energy, can hold 10 knots above the aircraft’s specified approaching speed; the flight path angle, which influences the potential part of the total energy, can be 1 dot higher (${\gamma }_{1dotup}$) than the standard gliding angle, which is 3$°$ for RWY 10 at Taipei Songshan Airport. As for the aircraft above the IMC gate, the aircraft’s deceleration performance is at 20 knots per nautical mile and still stays 1 dot above the 3$°$ glide path. The upper energy boundary is as follows:

$$E_{upper}=\left\{\begin{array}{c}D\times tan\left({\gamma }_{1dotup}\right)+\frac{{(V_{ref}+20)}^2}{2g} , when below IMC gate\\ D\times tan\left({\gamma }_{1dotup}\right)+\frac{{(V_{ref}+20+V_{cu})}^2}{2g} , when above IMC gate\end{array}\right.$$

(2)

where $V_{cu}=\left(D-D_{IMC}\right)\times 20\left(\frac{kts}{NM}\right)$ and D refers to the ground distance to the touchdown point.

In contrast to the upper boundary, the lower boundary is the minimum total energy of the aircraft when approaching. There is no difference in the definition of the lower boundary whether the aircraft is above or below the IMC gate. The potential part of the total energy, which is the aircraft’s gliding angle, is accepted at 1 dot lower than the glide path (${\gamma }_{1dotlow}$), while the speed can be 5 knots lower than the specified approaching speed. The lower energy boundary is as follows:

$$E_{lower}=\left\{\begin{array}{c}D\times tan\left({\gamma }_{1dotlow}\right)+\frac{{(V_{ref}-5)}^2}{2g} , when below IMC gate\\ D\times tan\left({\gamma }_{1dotlow}\right)+\frac{{(V_{ref}-5)}^2}{2g} , when above IMC gate\end{array}\right.$$

(3)

2.4.2. Ideal Energy State

The ideal energy state of an aircraft is the condition that the aircraft should have during the approach phase. When the aircraft is above the IMC gate, it should stay precisely on the 3-degree glide path, and the optimal descent rate is 10 knots per nautical mile. On the other hand, when the aircraft is below the IMC gate, the aircraft should be on the standard glide slope, $\gamma$, and the speed selected is the specified approaching speed ($V_{ref}$) according to the aircraft type, where $V_{ref}$ is derived from the manual of the corresponding aircraft type, which is the ATR 72-600 in this study. The ideal energy state is as follows:

$$\begin{gathered} E_{ideal}=\left\{\begin{array}{c}D\times tan\left(\gamma \right)+\frac{V_{ref}{}^2}{2g} , when below IMC gate\\ D\times tan\left(\gamma \right)+\frac{{(V_{ref}+V_c)}^2}{2g} , when above IMC gate\end{array}\right. \\ where V_c=\left(D-D_{IMC}\right)\times 10 \left(kts/NM\right) \end{gathered}$$

(4)

2.4.3. Total Energy State Observation

After the upper and lower energy boundaries and the ideal energy curve are defined (Figure 3), the aircraft’s total energy state during the approach phase can be observed, as shown in Figure 4, where the blue line is the upper boundary, the green line is the ideal energy curve, the orange line is the lower energy boundary, and the red line is the total energy of a specific flight. However, the actual total energy curve of the flight makes it difficult to observe the excess part compared to the energy boundary. Therefore, an energy normalization process must be applied to better study the flight to determine whether it is unstable or not [20]. The normalized total energy, the so-called energy factor (EF), can represent the deviation level of the aircraft’s energy state and is a function of the remaining distance to touchdown. The normalization process is shown in Equation (4), where $E_{dev}$ is the energy deviation between the actual flight energy and the ideal energy, $E_{d-margin}$ represents the difference between the ideal energy curve and the lower boundary, and $E_{u-margin}$ shows the difference between the upper boundary and the ideal energy curve. The normalized energy boundary is shown in Figure 5.

$$EF\left(D\right)=\left\{\begin{array}{c}\frac{E_{dev}}{E_{d-margin}} if E_{dev}<0\\ \frac{E_{dev}}{E_{u-margin}} if E_{dev}\ge 0\end{array}\right.$$

(5)

Figure 6 shows a flight’s normalized total energy on the normalized energy boundary. The blue line is the normalized upper boundary, the green line is the normalized ideal energy curve, the orange line is the normalized lower boundary, and the red line is the normalized total energy of a specific flight. After normalization, the tendency of the total energy can be observed easily, and a comparison of the deviation in the total energy with the boundary and ideal curve can be more straightforward.

Something that is noteworthy is that the selected window of the exceeded energy is different for energy deficit and energy excess. The exceeded total energy over the upper boundary is selected from the intermediate approach fix (IF) to the runway threshold while the total energy below the lower boundary is selected only after the flight passes the final approach fix since the approach speed limitation starts from six nautical miles away from the runway threshold, which is the outer marker (OM or FAF), and the pilot is mostly requested to maintain high airspeed down to the OM [4]. In Figure 7, the red dashed line shows the window of the selected area. As shown in Figure 8, the red area is the excess energy of the flight, which may cause an unstable approach and increase safety risks during landing. The excess area of the normalized total energy can be one of the indicators to identify whether the approach is unstable or not. However, the output that is put into model training is not the form of the excess area since the value of the excess area for each flight fluctuates too much and would reduce the model’s performance. The mean normalized total energy of the flight during the approach phase is the output instead, and the excess area would be an auxiliary to identify instability.

3. Data Preprocessing and Preparation

The ADS-B data downloaded from the FlightRadar24 website consist of data from the whole flight phase. As this study mainly considers the situation during the approach and landing phases, the flight data from the IF of Taipei Songshan Airport to the touchdown point were extracted, as shown in Figure 9. The corresponding METAR data were found and labeled by the ATMAP weather algorithm to classify the flight according to the weather conditions. On top of that, an energy boundary was built to define and identify stable approaches.

3.1. ADS-B Data Preprocessing

The original data downloaded from the website must be preprocessed before putting them into the model for training since some data are missing and incomplete. Outliers need to be eliminated so as not to affect the results.

3.1.1. Data Cleaning and Conversion

Since the ADS-B data originate from the FlightRadar24 website, the sampling rate of the data recorded is longer than the receiver that is self-installed, and sometimes, it will not record the data of the complete flight phase, so the record that lost the information during the approach and landing phases is cleaned. A total of 498 flights were pruned after data cleaning from the original 9422 flights. After that, the units or the representation of some parameters in the ADS-B dataset must be converted to make the calculation of the flight parameter inputs more convenient. The timestamp in the original dataset must be converted into UTC form so that the data’s time can easily be made clear. The vertical speed of the aircraft, which is in unit of feet per minute (fpm), needs to change to the unit of knots, which is a nautical mile per hour, for matching up with the ground speed when calculating the flight parameter inputs.

3.1.2. Landing Runway Identification

There are two runway numbers at Taipei Songshan Airport, namely runway 10 (RWY 10) and runway 28 (RWY 28), where the number represents the heading when the aircraft lands. Through the “direction” item in the ADS-B data of the landing flight, it can be separated into two landing directions, and landing runway identification can be carried out. Since the true directions of runway 10 and runway 28 are 92.22° and 272.24°, respectively, the flight, whereby the heading of the aircraft was maintained lower than 100 in the last few pieces of the ADS-B data, would be recognized as a flight landing on runway 10; otherwise, if the heading of the aircraft was maintained higher than 250, it would be recognized as a flight landing on runway 28. After data cleaning and landing runway identification, 8494 flights landed on RWY 10 and 1157 flights landed on RWY28 at Taipei Songshan Airport.

3.1.3. Outlier Detection

HDBSCAN is utilized for detecting abnormal flights that are spatial anomalies. The concept of a spatial anomaly is defined as a trajectory whose spatial metrics do not conform to an identified set of air traffic flows representing standard spatial operations [10]. HDBSCAN extends DBSCAN with the advantage of the ability to identify various density clusters since the flux of the traffic flow during the approach and landing phases might be different, i.e., some flows appear frequently while some are less frequent [10]. Since the outliers would be found typically after air traffic flow identification, the first step is to classify the trajectories of landing at Taipei Songshan Airport. Latitude and longitude are the two parameters used in HDBSCAN to identify the air traffic flow. As HDBSCAN is an unsupervised machine learning algorithm, it cannot define the genuine outliers; the proportion of outliers must be decided by modulating the HDBSCAN hyperparameter, the minimum cluster size. The heading of the landing flights after passing the altitude of 5000 feet, which is approximately 17 nautical miles radius from the touchdown point, is classified into three clusters, as shown in Figure 10, which are identified sets of air traffic flows representing standard spatial operations. The black line in Figure 10 represents the outliers where the minimum cluster size was set to 55 to detect five percent of outliers, which are 430 flights.

3.2. Weather Classification

The weather around an airport is an important factor affecting the air traffic, especially when aircraft take off or land. According to the statistics of the safety report 2021 edition produced by the IATA [5], the greatest contributing factor to accidents in the threat category is meteorology. However, this study wants to discover the unstable approach factor that is non-weather-related. As mentioned previously, the weather phenomena from the METAR of each flight were labeled by applying the ATMAP weather algorithm, as shown in

Table 3. Labeled weather phenomena.

Weather Phenomena	Visibility	Ceiling Height	Wind	Crosswind	Included Angle	Rain	Temp.	DP	Without	With
B78606 2019-09-27	4	3	3	1	1	2	3	3	3	4
B78722 2019-08-30	3	2	3	4	2	2	4	3	2	3

. Since visibility, ceiling height, and included angle are the top three weather conditions that affect RWY 10 at Taipei Songshan Airport the most [12], these three phenomena were the input features for the clustering task, which was performed by K-means.

The result of the weather phenomena clustering is shown in

Table 4. Clustering results by K-means.

	Cluster A			Cluster B			Cluster C			Cluster D
Flight Data Volumes	1267			5425			1241			561
TOP 3 Weight	Visibility	Ceiling Height	Included Angle	Visibility	Ceiling Height	Included Angle	Visibility	Ceiling Height	Included Angle	Visibility	Ceiling Height	Included Angle
Label	1	1	1	1	1	2	1	2	2	3	2	2
Mean Risk Index	8.3			9.8			10.7			12.6

. The silhouette score that was close to 0.7 was chosen, which is an acceptable result to explain the weather conditions around the airport, and the quantity of the cluster will not be too much. The flights were divided into four clusters, with the number of flights being 1267, 5425, 1241, and 561. In cluster A, the labels for visibility, ceiling height, and included angle are all 1, which is the cluster consisting of the flights with the lowest weather risk index. For cluster B, only the included angle is labeled with 2, while the others are 1. As for cluster C, ceiling height and included angle are both tagged by 2, and visibility is still 1. Finally, the cluster with the highest weather risk index is cluster D, with visibility labeled with 3 and both ceiling height and included angle labeled with 2. After weather classification, there were four DL models built according to different weather conditions to analyze the aircraft’s approach risk factors in the same weather situation.

3.3. Flight Parameters

When building a deep neural network model, several features of the aircraft are required as the model input and output to train the model. In [21,22], the features used to train the random forest model were derived from the FDR. Since there are more than six hundred parameters in the FDR, parameter selection must be carried out first to ensure the parameters used can advance the model’s training more efficiently and faster. The features used to train the model in this paper can be obtained from nothing but the ADS-B data sources, which limited the selection of the parameters since there are few basic parameters in ADS-B. The parameters were chosen according to the results and discussions in [21,23], where the input parameters were categorized into three groups: aircraft historical performance parameters, energy-related parameters, and trajectory-related parameters, as shown in

Table 5. Input flight parameters.

Categories	Item
Historical Performance	Altitude at level flight
	Time spent in level flight
	Speed at level flight
Energy-Related	Speed	at 3000 ft = 1000 ft + 1000 fpm × 2 min
	Rate of descent
	Horizontal distance
	Air miles to touchdown
Trajectory Related	Lateral deviation
	Vertical deviation
	Flight path angle deviation

.

3.3.1. Aircraft Historical Performance Parameters

The features in this category represent the aircraft’s state during the level flight phase. These historical features include altitude at level flight, time spent in level flight, and the speed at level flight. Since level flight is a long-lasting flight phase, the parameters derived were the mean value of the whole phase.

3.3.2. Energy-Related Parameters

The energy-related parameters represent the aircraft’s energy state, including kinetic and potential energy states. While horizontal distance and air miles to touchdown may not directly seem energy-related, their relationship with specific altitude correlates them to the potential energy state. To cite an instance, a shorter horizontal distance to the touchdown point at a higher altitude indicates a higher energy state of the aircraft in the final approach phase, which means a more considerable possibility of an unstable approach happening. The parameters of the aircraft’s kinetic energy state enclose the ground speed and the rate of descent of the aircraft at a specific altitude, which are correlated with the speed and can be obtained from the ADS-B data directly.

3.3.3. Trajectory-Related Parameters

Lateral deviation, vertical deviation, and flight path angle at a particular altitude are the aircraft’s trajectory-related parameters for the inputs of the NN model training. Lateral deviation and vertical deviation are the lateral and vertical distance between the aircraft’s position and the standard glideslope. The flight path angle is calculated from the ground and vertical speeds, as shown in Equation (6). The trajectory deviation, the distance between the aircraft position and the 3$°$ glide slope, of the flight is one of the model outputs.

$$Flight Path Angle={\mathit{tan}}^{-1}\left(\frac{Vertical Speed \left(VS\right)}{Ground Speed \left(GS\right)}\right)$$

(6)

After the model has been generated, these input flight parameters of the flights in one of the weather clusters will be trained for learning to predict the outputs, which are the trajectory deviation, and another one will be discussed later.

4. Model Training and Discussion

The deep neural network model used for predicting the unstable approach is introduced in this section, including the description of the data used for model training, the model architecture, the training results, and the discussion of the non-weather-related factors analysis.

4.1. Overview of the Deep Neural Network

4.1.1. Input/Output Data Selection

Table 6. Input and output data of the model.

	Item	Attribute
Input	Altitude at level flight	Value
	Time spent in level flight	Value
	Speed at level flight	Value
	Speed	Value
	Rate of descent	Value
	Horizontal distance	Value
	Air miles to touchdown	Value
	Lateral deviation	Value
	Vertical deviation	Value
	Flight path angle	Value
Output	Normalized total energy	Value
	Mean trajectory deviation	Value

shows the data involved in model training, including the input and output data. There are ten features used for model training input, and all of them are in numerical values. The model training outputs are the flight’s normalized total energy and the mean trajectory deviation during the approach and landing phases, which starts from the IF to the touchdown point.

4.1.2. Data Interpolation and Normalization

Interpolation is a technique to estimate and create new data points from a known dataset for discrete functions in numerical analysis when there is a limited amount of data. Unlike data fitting, the new curve with new data points created by interpolation should pass through the original data. Owing to the sparsity of the ADS-B data downloaded from the FR24 website, there are different numbers of data points for each flight during the phase that we are interested in. Since the calculation of the difference between the normalized total energy and the normalized energy boundary, which is the excess area, is an integral process, more data points can obtain more precise results. Figure 11a,b shows the flight data before and after the interpolation process.

Due to the various range and units of the input parameter values, it is essential to rescale these parameters to a united dimensionless form for optimizing the model and increasing the precision of the prediction. Data normalization is a technique to scale the original data into intervals of 0 to 1 and not change the data’s original distribution. There are some common normalization methods, such as standard deviation normalization and min–max normalization, where the type used in this paper is the latter, as shown in Equation (7), where $X_{max}$ and $X_{min}$ are the original maximum and minimum value in that parameter, respectively, and $X$ and $X_{norm}$ are the value before and after normalization, respectively.

$$X_{norm}=\frac{X-X_{min}}{X_{max}-X_{min}} \in \left[0,1\right]$$

(7)

4.1.3. Hyperparameter Tuning

In each machine learning algorithm, some parameters that need to be settled on before or after the learning process are called hyperparameters. These hyperparameters can control the performance of the model learning process, such as the learning efficiency and the prediction accuracy, etc. Take the machine learning model used in this paper as an example, the type of activation functions, the number of hidden layers, the number of neurons in each layer, etc., are the hyperparameters of the NN model, and all must be decided in advance.

Since there might be lots of combinations of hyperparameters, a method called grid search with cross-validation (GridSearchCV) was applied to help select the best combination of the hyperparameters to enhance the model’s performance. GridSearchCV is a hyperparameter tuning method that is composed of grid searching and cross-validating, where grid search means using the method of exhaustion in a series of parameters and cross-validation means finding the candidate that makes the most precise results. This study makes use of GridSearchCV to tune and determine the suitable number of neurons for the input layer and two hidden layers in the NN model. The result of the hyperparameter tuning process, which is measured by mean squared error (MSE), for clusters A~D is shown in

Table 7. Results of GridSearchCV.

	Cluster A	Cluster B	Cluster C	Cluster D	All Flights
MSE Score	0.0075	0.0039	0.019	0.0068	0.0036
Number of Neurons for Each Layer	(22, 23, 6)	(22, 17, 7)	(11, 12, 7)	(20, 21, 10)	(20, 21, 6)

. The mean squared error shows the difference between the actual value and the predicted value, while the combination with the lowest MSE score would be the one used for model training.

4.1.4. Model Architecture

The model used for the prediction in this study was constructed by a deep neural network, which is composed of an input layer with several input parameters, some hidden layers, and an output layer with output features as a function of the input features. The outputs of the model can be either numerical type for regression NN or categorical type for classification NN. In this paper, the goal was to predict the normalized total energy and the trajectory deviation of the aircraft, which are numerical predictions, and thus the regression model is to be implemented.

The training model was designed for 10 flight parameter inputs and 2 outputs (mean normalized energy and trajectory deviation), which can explain the aircraft’s energy state and the potential for an approach to become unstable. There is one input layer, two hidden layers, and an output layer for each model. These layers in the neural network model are fully connected, and the input parameters will pass through these layers. Then, these layers follow an activation function to activate the outputs of each layer, where the rectified linear unit (ReLU) is the activation function used in the hidden layer in this paper. After each layer, a batch-normalization process is implemented, which can reset the distribution of the scattered data, help solve the gradient vanishing problem, and speed up the result to convergence.

The architecture of the cluster A model is shown in Figure 12, where the numbers of neurons for the four layers in cluster A are 22, 23, 6, and 2. With the same architecture of the cluster A model, the numbers of neurons for the four layers in cluster B are 22, 17, 7, and 2; the numbers of neurons for the four layers in cluster C are 11, 12, 7, and 2, and the numbers of neurons for the four layers in cluster D are 20, 21, 10, and 2.

4.2. Model Training and Testing Results

The training data used in the proposed model mentioned earlier are the ADS-B data of ATR 72-600 flights landing on RWY 10 at Taipei Songshan Airport. There are approximately 9000 flights in a year (2019/07~2020/06), which are split into four groups (clusters A~D) in the weather classification step. Furthermore, the data are split up into 80 percent of the training dataset and 20 percent of the testing dataset, as shown in

Table 8. Training and testing data volume of each cluster.

	Cluster A	Cluster B	Cluster C	Cluster D	All Flights
Training Data	1014	4340	993	449	6795
Testing Data	253	1085	248	112	1699

. After the NN model is trained, there is a loss function, which is the MSE in this model, evaluating the model performance. MSE considered the difference between the predicted value estimated by the model and the true value from the testing dataset.

The NN model’s forecast accuracy, which is a measurement of the predicting model’s performance, is calculated by using the weighted average percentage error (WAPE) in the terms of percentages. The WAPE, or MAD/mean ratio, is a common validation method for regression models, such as the NN model in this paper. It is the sum of the absolute error divided by the sum of the actual values and therefore gives more importance to errors, as shown in Equation (8) where $A_t$ is the actual value, $F_t$ is the forecast value, and n represents the total data number. The difference between the WAPE and the MAPE, which is another commonly used validation measure, is that the WAPE can counter the situation when there are zeros or values close to zero in the denominator. As the WAPE is the error of the model, the smaller it is, the better. Afterward, the accuracy of the model prediction in Equation (9) is defined as one hundred percent minus the WAPE as the result.

$$WAPE=\frac{{{\displaystyle \sum }}_{t=1}^n\left|A_t-F_t\right|}{{{\displaystyle \sum }}_{t=1}^n\left|A_t\right|}\times 100\%$$

(8)

$$Accuracy=100\%-WAPE$$

(9)

The training results of cluster A, which contains the smallest number of weather risk index flights, are displayed in Figure 13, Figure 14 and Figure 15. The test loss of the model for cluster A is approximately 0.01, and the accuracy of the prediction for the normalized total energy and trajectory deviation is 87.9% and 81.8%, respectively, which means the WAPE is 12% and 18% respectively, which is an acceptable result.

Cluster B and cluster C are the two clusters with middle-weather risk index flights, where cluster B comprises most of the flights landing on RWY 10. The test losses of cluster B and cluster C are 0.00364 and 0.0172, respectively. The accuracy of the prediction for the normalized total energy and trajectory deviation is 84.16% and 81.19% for cluster B and 83.16% and 84% for cluster C, respectively.

The results of the model consisting of all the flights landing at RCSS RWY 10 are listed in

Table 9. Model training and testing results.

	Cluster A	Cluster B	Cluster C	Cluster D	All Flights
Training Loss	0.00834	0.00333	0.01563	0.00288	0.00317
Test Loss	0.00822	0.00364	0.01723	0.00595	0.00310
Normalized Total Energy Accuracy	87.93%	84.16%	83.15%	87.19%	85.15%
Trajectory Deviation Accuracy	81.80%	81.19%	84.01%	81.68%	82.11%

. The test loss of the model for all flights is approximately 0.003, and the accuracy of the prediction for the normalized total energy and trajectory deviation are 85.15% and 82.11%, respectively. These results have shown that in different kinds of weather conditions, or not considering the weather, the models have similar good performance.

4.3. Unstable Approach Identification

After the value of the normalized total energy and the trajectory deviation are predicted, they can be used for identifying whether the approach of that flight is unstable. If either of these two values exceeds the threshold set in advance, the approach would be recognized as unstable. The thresholds of identifying an unstable approach by normalized total energy and trajectory deviation depend on the median and quartile of each cluster, as shown in

Table 10. Unstable detection threshold for each cluster.

	Cluster A	Cluster B	Cluster C	Cluster D	All Flights
Normalized Total Energy Threshold	1.56	1.76	1.61	1.75	1.715
Trajectory Deviation Threshold (ft)	980	985	1085	1165	1010

. The excess area of the energy curve would be a confirmation. When the excess area exceeds a specified value, the approach has a higher possibility to be identified as an unstable situation.

In Figure 16 and Figure 17, the boxplot shows the normalized energy and trajectory deviation distribution of all the flights landing on RWY 10 at RCSS. The orange line shows the median of the whole dataset and both sides of the box are the lower quartile ($Q_1$) and the upper quartile ($Q_3$), respectively. The small circles in the boxplot are outliers where their distance from $Q_1$ or $Q_3$ is more than $1.5\times IQR$ ($\mathrm{where}IQR=Q_3-Q_1$). Therefore, the threshold of normalized total energy and trajectory deviation is 1.715 and 1010 feet, respectively, to detect unstable approaching flights.

Table 11. Outlier distribution.

	Cluster A	Cluster B	Cluster C	Cluster D	All Flights
Number of Flights	1267	5425	1241	561	8494
Outliers	83 (7.6%)	436 (8.2%)	121 (8.7%)	72 (8.4%)	712(8.4%)

shows the number and percentage of the outliers in each cluster by applying the threshold obtained from all flights, where the outlier proportions in each cluster are all below 10%.

4.4. Feature Importance Analysis

Feature importance analysis is a crucial step in identifying the importance of different parameters affecting approach safety. As mentioned earlier, feature importance is a measurement to present the importance of the input parameters depending on the usefulness when making the prediction. In this work, permutation feature importance was selected to calculate the input parameter weights of the model for the four clusters. It can be seen in

Table 12. Permutation feature importance for clusters A~D.

	Cluster A	Cluster B	Cluster C	Cluster D
Altitude at level flight	0.04634	0.04923	0.06502	0.02321
Time spent in level flight	0.06249	0.06289	0.07309	0.02796
Speed at level flight	0.04787	0.05156	0.07241	0.04565
Speed at 3000 ft	0.08524	0.10213	0.11418	0.05124
Rate of descent at 3000 ft	0.03405	0.03262	0.04830	0.02260
Horizontal distance at 3000 ft	0.08887	0.13008	0.07581	0.06856
Air miles to touchdown at 3000 ft	0.453308	0.40496	0.35348	0.41451
Lateral deviation at 3000 ft	0.05618	0.05833	0.06034	0.02483
Vertical deviation at 3000 ft	0.07152	0.05695	0.06812	0.07715
Flight path angle at 3000 ft	0.05436	0.05124	0.06923	0.04029

that the feature, air miles to touchdown point when the aircraft is at 3000 ft, is the most important feature of the models for all air traffic flows. The speed and the horizontal distance at the IF have a moderate effect in clusters A, B, and C. On the other hand, the model results of cluster D are affected by the vertical deviation and horizontal distance more at the IF. As stated in [21], the features that have higher feature importance at a higher altitude can likely be determined as the precursors to the unstable approach event, which could be an alert in advance. The air-mile to touchdown indicates that the potential energy state is of the greatest importance for flights landing at RCSS RWY 10. The speed at 3000 ft ranks second for the feature in clusters A, B, and C and correlates to kinetic energy. The results show that energy management is crucial for approach safety since the top three important features for most flights are energy-related factors.

4.5. Comparison of Approach Risk Analysis with Weather Data and Energy Management

A previous study [12] showed the evaluation of approach risk considering weather conditions. Although weather conditions have a significant influence on approach phase safety, there are still the other non-weather-related factors, such as pilot operations. To study the effects of non-weather factors, energy management was introduced in this study to evaluate the energy factors of flights with similar risk analysis results and weather conditions. In [12], the authors considered weather as one of the important factors affecting the safety of the approach phase and gave a risk index to each flight by using the weather labels and the corresponding weights.

In this study, the flights were clustered into four groups according to the risk index so the flights in the same group are assumed to have similar weather conditions during the approach phase. If the flights in the same cluster have different results of unstable approach analysis, the flight that is determined as unstable has a greater possibility to be affected by non-weather-related factors. Taking the flights in cluster A as an example, cluster A is the cluster that has the lowest mean weather risk index compared to the other clusters. There are some flights identified as unstable since they have higher specific total energy and have an excess area of the energy curve exceeding the threshold. In Figure 18, flight AE1276 on 3 September 2019 has a weather risk index of 8.97, which is in the middle-risk range but is identified as unstable by the proposed energy management method. On the other hand, flight AE366 on the same day has a higher weather risk index when approaching which is 9.49, but it is identified as a normal landing flight through the energy observation as shown in Figure 19. Therefore, for the results, there could be some non-weather-related factors that affected flight AE1276 when approaching and landing.

There are other flights that have similar results, such as flight AE374 on 19 May 2020, which has a low-risk weather index of 7.24. However, it was identified as an unstable flight regarding energy management since the mean normalized total energy was up to 1.1. As shown in Figure 20, the energy excess part is close to the minimum stabilization height (1000 ft). If the go-around is required for this flight, the risk will increase over time since it needs more time to recover the altitude and the airspeed of the aircraft during go-around.

5. Conclusions

This study proposed a flight safety analysis approach during the approach and landing phases regarding energy management to figure out the factors that could affect flight safety. The energy calculated from the ADS-B data was applied to determine unstable approaches and to find out the non-weather-related factors that affect the flight safety of approach. To evaluate the influence of weather- and non-weather-related factors, the outliers of the flights classified by weather labels were detected and eliminated from the analysis by applying HDBSCAN. The deep learning method was utilized to predict the normalized total energy and trajectory deviation of the flights at Taipei Songshan Airport. The prediction accuracy for the normalized total energy and trajectory deviation of all flights were higher than 82.15%. The results show that in different kinds of weather conditions, or not considering the weather, the models have similar good performance. From the results, there are four clusters which denote the different weather situations of the flights. Finally, the discussions and comparisons about the non-weather-related factors and unstable approach flights were conducted to discuss the influence of non-weather-related factors that make the approach unstable. From the results, it can be seen that the energy management of the aircraft provides different points of view for flight safety analysis for the flights during the approach phase.

Despite the ADS-B data being incomplete for human factor safety analysis, this study provides the feasibility to analyze the factors that affect approach safety through energy management. As pointed out previously, it would be better if the ADS-B data were received by the owned device since the ADS-B data downloaded have a disunity quantity or data missing problem which leads to a trade-off in data volume. To find out more specific causes of unstable approaches, the data in the flight data recorder will be required in model training in the future, and then it will be possible to go one step further to increase the detection accuracy of unstable approaches and to discover the pilot errors related to unstable events. Furthermore, the model built in this study could be applied to other airports that have the same landing aircraft category as the ATR 72-600 since the concept of analyzing unstable approaches with energy management is not limited by the airport conditions or airline policies. Finally, a validation process is essential to prove the model’s practicality. Nonetheless, there are no data of recent approach accidents or events that can be used for model validation in this study. It will be a must-do in the future to collect the real event data to validate and to improve the model, which will make the post-processing analysis more comprehensive.