1.2.1 VIL

VIL is used to describe the amount of liquid water accumulated in the entire atmospheric column from the altitude of the radar to the top of the atmosphere. It is calculated as shown in Equation (2). Higher VIL values usually indicate the presence of a significant amount of liquid water within a thunderstorm, often foreshadowing the occurrence of severe convective weather such as heavy rainfall, hail and strong winds.

(2) \begin{align} {\rm{VIL}} = {\rm{\;}}\sum 3.44 \times {10^{ - 6}}{\left( {\frac{{{Z_i} + {Z_{i + 1}}}}{2}} \right)^{4/7}}\Delta h, \end{align}

where ${Z_i}$ is the radar reflectivity and Δh is the height expressed in meters.

Table 1.

Convective weather classification

In this paper, the VIL data are stored in mat files, and each element in the matrix represents the VIL value of 0.01°^*0.01° airspace cell at the specified location. In order to express the severity of convective weather in the airspace, this paper has counted the VIL percentile value of the airspace for the specified segment. Furthermore, in order to measure the coverage of convective weather in the airspace, this paper has counted the proportion of the number of VIL ≥ 3.5 kg/m² rasters in the airspace of the specified segment. Therefore the VIL features used in this paper are:

1. 1. 50 ^th VIL: the 50th percentile VIL value in the airspace of the segment;

2. 2. 75 ^th VIL : the 75th percentile VIL value in the airspace of the segment;

3. 3. 90 ^th VIL : the 90th percentile VIL value in the airspace of the segment;

4. 4. maxVIL : the maximum VIL value in the airspace of the segment;

5. 5. Proportion of VIL ≥ 3.5 kg/m²: the number of grids with VIL ≥ 3.5 kg/m² as a percentage of all grids in the airspace of the segment.

1.2.2 ET

ET is the height corresponding to an echo of 18.3 dBZ detected by meteorological radar, as shown in Equation (3). It is used to measure the intensity and vertical development of convective activity. Higher ET values usually imply higher convective cloud tops, indicating that severe convective weather might be developing or has already developed. The ET features used in this paper include:

1. 6. 50 ^th ET: the 50th percentile ET value in the airspace of the segment;

2. 7. 75 ^th ET: the 75th percentile ET value in the airspace of the segment;

3. 8. 90 ^th ET: the 90th percentile ET value in the airspace of the segment;

4. 9. maxET: the maximum ET value in the airspace of the segment.

(3) \begin{align}{\rm{ET}} = {{\rm{h}}_{BR \gt threshold}} \end{align}

1.2.3 CR

The basic reflectivity (BR) after logarithmic transformation of Z can reflect the magnitude of Z value as in Equation (4).

(4) \begin{align} BR = 10{\rm{log}}\left( Z \right) \end{align}

CR refers to a product derived within a single volume scan that projects the maximum reflectivity factor found in constant elevation azimuth scans onto Cartesian grid points. It represents the maximum reflectivity factor within a unit volume, as depicted in Equation (5).

(5) \begin{align} CR = maxBR \end{align}

A larger value of CR represents a higher convective intensity.

1. 10. 50 ^th CR: the 50th percentile CR value in the airspace of the segment;

2. 11. 75 ^th CR: the 75th percentile CR value in the airspace of the segment;

3. 12. 90 ^th CR: the 90th percentile CR value in the airspace of the segment;

4. 13. maxCR: the maximum CR value in the airspace of the segment;

5. 14. Proportion of CR ≥ 41 dBZ: the number of grids with CR≥41dBZ as a percentage of all grids in the airspace of the segment.

1.2.4 Longitudinal weather features of segment

This article analyses longitudinal and lateral convective weather conditions in sector route segment to identify their impact on air traffic flow.

The impact of longitudinal weather in sector route segment on air traffic flow is reflected in the length of the path of a flight crossing the convective weather portion, indicating the severity of longitudinal convective weather in sector route segment. Figure 3 presents the method for assessing longitudinal convective weather in sector route segment.

Figure 3.

Schematic of longitudinal weather feature.

Parallel lines are drawn along the direction of the sector route, spaced at intervals of 1 km parallel to the centerline. The VIL/CR values of all grids along each parallel line are calculated. When the VIL/CR value exceeds level 3 weather, the raster is considered to be in convective weather, and the convective weather weight value of the raster is set to be 1; otherwise, the weight value of the raster is VIL/3.5 or CR/40. The sum of the convective weather weights of all the rasters on each parallel line indicates the longitudinal weather severity of the line as shown in Equation (6).

(6) \begin{align}lin{e_j} = \mathop \sum \limits_{i = 1}^n {w_{ij}},{\rm{\;\;}}{w_{ij}} = \left\{ {\begin{array}{*{20}{c}}{1,\;\;if\;VIL \geqslant 3.5\;\;or\;CR \geqslant 40}\\[5pt] {\frac{{value}}{{3.5}}or\frac{{value}}{{40}},\;\;if\;VIL \lt 3.5\;\;or\;CR \lt 40}\end{array}} \right.,\end{align}

Where $\;{w_{ij}}$ denotes the convective weather weight of the i ^th grid of the j ^th parallel line on the segment. $Lin{e_j}$ denotes the convective weather duration of the flight on the j ^th parallel line, and value denotes the VIL or CR value of each grid on the parallel line in this function.

The longitudinal route segment weather features are:

15/16. VIL_parallel_mean/CR_parallel_mean: The mean value of the impact along each parallel line represents the longitudinal average impact (based on VIL/based on CR);

17/18. VIL_parallel_max/CR_parallel_max: The maximum value of impact among all parallel lines (based on VIL/based on CR);

19/20. VIL_parallel_min/CR_parallel_min: The minimum value of impact among all parallel lines (based on VIL/based on CR);

21/22. VIL_parallel_median/CR_parallel_median: The median value of impact among all parallel lines (based on VIL/based on CR).

1.2.5 Lateral weather features of segment

When convective weather does not entirely cover the cross-section of a sector route segment, flights along the route can still traverse areas not affected by the weather. The lateral impact of convective weather on the sector route segment manifests in the degree of obstruction in the direction of the route’s width, quantifying the proportion of the route’s width covered by convective weather in the vertical direction, as indicated in Equation (7).

(7) \begin{align} {W_r} = \frac{{mcw}}{{width}} \end{align}

where mcw denotes the minimum available width on the route segment, width denotes the route segment width, and Wr represents the available proportion of traffic flow within the sector route segment.

The determination of the minimum available width of the segment, as illustrated in Figures 2, 3 and 4, begins by identifying the grid cells within the sector route affected by convective weather; Subsequently, the DBSCAN algorithm is employed to cluster multiple convective weather blocks within the route, depicting the progression from Fig. 4(a) and (b); The clustered convective weather blocks are utilised with a convex hull algorithm to generate a convex hull region, denoting the area impassable for aircraft, depicting the process from Fig. 4(b) and (c).

Figure 4.

Flow percentage calculation process diagram for crossable sector route segments.

Finally, the Dijkstra algorithm is applied to determine the mcw in the vertical direction within the sector route. The lateral features of the flight path in this section are included:

1. 23. Cross_flow_VIL: the proportion of traffic flow that can traverse a route, calculated based on VIL;

2. 24. Cross_flow_CR: the proportion of traffic flow that can traverse a route, calculated based on CR.

3.2.1 Label

This paper uses airspace availability as the label, airspace availability denoted as SRA _r, which is calculated as shown in Equation (8).

(8) \begin{align} SR{A_r} = \frac{{SR{C_r}}}{C} \end{align}

In Equation (8), SRC _r and C represents the actual flight number under convective weather and the maximum flight number under clear weather of the sector route in the corresponding period. The number of flights in sector route segment varies significantly across different time intervals, typically with fewer flights in the early hours and a higher volume during daylight hours. However, there exists a degree of similarity in the number of flights during equivalent time segments on different days. This study involves the calculation of sector route segment capacities for each hour within a day over a specific time span (usually a month), as derived in Equation (9).

(9) \begin{align} {C^p} = {90^{{\rm{th}}}}\left\{ {C_1^p,C_2^p, \ldots C_i^p} \right\} \end{align}

where ${C^p}$ represents the hourly capacity of a sector route segment during a specific period of a day, and $C_i^p$ represents signifies the sector route segment’s hourly capacity for a particular day within the statistical time span. Within a given time span, there might be instances where a day’s sector route segment hourly capacity exceeds its limit. We take the 90th percentile segment capacity value to represent the maximum accommodatable flight volume for that specific hour within the sector route segment. The availability label is the current hourly flight traffic divided by the current hourly capacity.

In the later part, SRA _p and SRC _p are the predicted value of availability and capacity, respectively.

3.2.2 Model evaluation metrics

In the field of machine learning, a common metric for assessing algorithm accuracy is the generalisation error. However, directly calculating the generalisation error is not practical in applications due to the unknown distribution of the sample data. Therefore, we utilise a validation set constructed from sample data and compute various performance metrics to estimate the generalisation performance of the RF algorithm.

This study employs four commonly used regression prediction metrics, including mean squared error (MSE), root mean square error (RMSE), mean absolute error (MAE) and R². The calculations are depicted in Equations (10), (11), (12) and (13).

(10) \begin{align} MSE = {\rm{\;}}\frac{{\mathop \sum \nolimits_{i = 1}^n {{\left( {{y_i} - {f_i}} \right)}^2}}}{n} \end{align}

(11) \begin{align} MAE = {\rm{\;}}\frac{{\mathop \sum \nolimits_{i = 1}^n \left| {{y_i} - {f_i}} \right|}}{n} \end{align}

(12) \begin{align} RMSE = \sqrt {\frac{{\mathop \sum \nolimits_{i = 1}^n {{\left( {{y_i} - {f_i}} \right)}^2}}}{n}}, \end{align}

(13) \begin{align} {R^2} = 1 - \frac{{\sum {{\left( {{y_i} - {f_i}} \right)}^2}}}{{\sum {{\left( {{y_i} - \bar y} \right)}^2}}} \end{align}

MSE measures the mean of the squares of the prediction errors. A smaller MSE value indicates that the model’s predictions are closer to the actual values, signifying better model performance. RMSE is the square root of MSE. It provides the standard deviation of the prediction errors. RMSE offers an intuitive feel for the size of the errors; compared to MSE, RMSE has the same units as the original data, making it easier to understand. MAE provides the average size of the prediction errors and, unlike MSE and RMSE, does not assign higher weight to larger errors. R² is a metric that measures the model’s ability to explain variability in the data. An R² value closer to 1 indicates that the model explains a higher variability, leading to better predictive performance.

3.2.3 Model training methods

This study utilises K-fold cross-validation method to assess the model’s generalisation performance and subsequently adjust the RF parameters. Cross-validation can effectively mitigates the risk of overfitting and accurately reflects the model’s generalisation capabilities. The RF design in the CWSRC model is illustrated in Algorithm 1. Within the K-fold cross-validation method, the original dataset is randomly partitioned into k equally sized subsets. Each iteration selects one subset the model’s validation dataset while using the remaining k−1 subsets as the training data. This process is repeated k times, rotating through each subset as the validation set, until all subsets have been used for validation. The metrics results obtained from each training validation are recorded, and the average of these k results yields the overall training performance for a particular set of hyperparameters. This iterative training and testing process makes the model performance assessment more reliable and reduces the sensitivity to the data partitioning method.

Algorithm 1.

RF.

3.2.4 Feature selection

This investigation incorporates 24 comprehensive features, which, while enriching the dataset, also extend the data collection phase and augment the model’s complexity. Effective feature selection is pivotal, as it not only simplifies the model by eliminating features with high intercorrelations but also enhances its generalisation capabilities and abbreviates the training duration. Conventionally, RFE operates by sequentially removing the least significant feature. This procedure, although systematic, may inadvertently undermine the model’s performance, particularly when the excised feature plays a pivotal role in certain combinations.

To address these shortcomings, the present study proposes an advanced modification of RFE that integrates aspects of stepwise regression. Following the exclusion of each feature, a comprehensive assessment is conducted to evaluate its interaction with other features and overall impact on model performance. Should this assessment reveal that certain eliminated features substantially bolster the model’s efficacy, these features are subsequently reinstated in the feature set during later iterations. This meticulous approach refines the feature selection process, culminating in an optimal feature subset that significantly enhances the model’s accuracy and robustness. Algorithm 2 outlines the algorithmic process.

Algorithm 2.

Recursive Feature Elimination Algorithm

In Algorithm 2, further stepwise explain are as follows:

1. (1) Initialisation (Lines 1–3): Initialise ‘Indicator’ to store performance indicators and ‘current_number_of_features’ to track feature count. Initialise the ‘recursiveFeatureElimination’ function.

2. (2) Loop and Training (Lines 4–5): Enter a loop to iteratively eliminate features until ‘current_number_of_features’ reaches ‘n_features_to_select’. Train the model (‘model.fit’) on the current feature set (‘x’ and ‘y’).

3. (3) Feature Importance Calculation (Line 6): Calculate feature importance using the model’s ‘feature_importances_’ attribute.

4. (4) Feature Removal (Lines 7–10): Identify and remove the least important feature from ‘x’ based on calculated importance. Update ‘current_number_of_features’ accordingly.

5. (5) Performance Evaluation (Lines 11–13): Evaluate model performance using the reduced feature set (‘x’ and ‘y’). Record the performance indicator (‘indicator’) in ‘Indicator’.

6. (6) Comparison and Recording (Lines 14–17): Compare the current performance indicator with the last recorded one (‘Indicator[−2]’). If performance worsens, record the removed feature as ‘delete_feature’.

7. (7) Recursive Call (Lines 18–19): Recursively call ‘recursiveFeatureElimination’ to continue feature elimination with the updated feature set.

8. (8) Return (Line 20): Once ‘current_number_of_features’ equals ‘n_features_to_select’, return the subset of features selected (‘subset = [current_feature, delete_feature]’).

Begin by establishing the desired number of target features. the original features is then input into the RF algorithm, the results of the performance metrics obtained from this training will be recorded, and the importance of the features will be ranked according to the RF’s feature importance evaluation function. Afterwards, selecting features based on the importance derived from random forest can identify the most influential features on the model’s predictive outcomes. After eliminating the features with the least importance, proceed to the next round of training with the remaining features. Record any excluded features when the metrics from the subsequent training round show a decline compared to the previous results. The process is repeated until the number of features reaches the target number of features, and the recorded excluded features are added to the remaining features to obtain the selected subset of features after screening.

3.2.5 Hyperparameter determination

The number of decision trees is a crucial hyperparameter in the RF algorithm. RF enhances the generalisation and stability of the model by aggregating the predictions of multiple decision trees. Each decision tree conducts random sampling on the training data and performs splits on randomly selected subsets of features, mitigating overfitting risks. Thus, while increasing the number of decision trees in a RF typically improves performance, there’s a limit to the increase as adding trees indefinitely may not guarantee improved performance. Considering the cost of computational resources and time consumption, it’s essential to establish a range of decision tree numbers. In this research, we follow precedents set in prior studies by setting the number of trees in the RF algorithm between 0 and 500 [Reference Genuer, Poggi and Tuleau-Malot19]. By traversing this range of parameters and observing the relationship between performance metrics and the number of decision trees through visualisations, we can determine the trend of performance metrics’ variations and ascertain the optimal number of decision trees.