Machine learning evaluation model of pilot workload in a low-visibility environment

ECG data preprocessing

ECG signals are weak, with an amplitude of 10 uV–4 mV and a frequency of 0.05–100 Hz. Moreover, ECG signals are easily influenced by external interference and pilots’ own actions, resulting in power frequency and electromyographical interference²⁷. In this study, ECG signals were preprocessed by wavelet transform, which could decompose the original signal into high- and low-frequency coefficients that were then denoised and reconstructed²⁸. The specific process is as follows:

The signals were processed via discrete wavelet transform, namely,

α and β in the above formula were discretized, then \(\alpha ={2^{ – j}},\beta ={2^{ – j}}k\left( {j,k \in Z} \right)\).

On this basis, discrete wavelet transform could be acquired as follows:

where

In this study, the ECG signals were decomposed using the Haar wavelet, as defined below³⁰:

where the scale function is

Wavelet denoising was performed through the soft threshold method³¹, i.e.,

where \({\omega _{j,k}}\) is the wavelet decomposition coefficient, and λ is the selected threshold.

Through the above method, filtering was performed via MATLAB, and the waveforms before and after processing are exhibited in Fig. 2. The processed ECG signals were subjected to QRS wave positioning, the R-wave position was marked, and the difference between two adjacent R waves, i.e., RR interval, was solved, followed by time- and frequency-domain nonlinear analysis.

Then, ECG R waves were detected using the Pan–Tompkins algorithm²⁹. First, the set threshold was compared with the maximum value of the signal. If the latter was greater than the former, the maximum value was defined as R-wave crest, and the threshold was timely updated. After the signal threshold was continuously updated, the specific position of ECG R wave could be accurately detected, effectively reducing the false detection rate. The detection method is as follows³⁰.

1. (1)

   The maximum value of the preprocessed ECG signal was extracted as the reference point of R-wave crest, and the signal threshold ts and noise threshold tn were initialized through the following formulas:

where pEEG is the preprocessed ECG signal, and fs represents the sampling frequency.

1. (2)

   The R-wave crest sample was compared with the signal threshold. If the former was greater than the latter, the position of this point was recorded, and the signal level sl was updated; if the signal threshold was greater than the R-wave crest sample, this point would be taken as a noise point, and the noise level nl was updated, specifically as follows:

where mp is the R-wave crest sample.

1. (3)

   The signal threshold and noise threshold were updated, as shown below:

The above steps were repeated, and the detected R-wave crest position was recorded. When the R-wave crest was detected using the above formulas, if the interval between the current and previous R-wave crest positions exceeded 1.66 times of the interval of the first eight R-wave crests, then missed detection occurred. In this case, the signal and noise thresholds were reduced to the half of the original thresholds, and redetection was performed from the previous R-wave crest position. In the process of backtracking detection, the signal and noise thresholds are expressed as

After preprocessing the ECG signal, R-wave peaks were identified using a combination of amplitude thresholding and peak detection algorithms. Specifically, a peak-prominence threshold and a minimum inter-peak interval of 600 ms were applied to avoid physiologically implausible detections. Peaks with amplitudes below the mean signal level plus 0.3 units were excluded during post-processing to reduce the risk of spurious detections.

As an illustrative case, the filtered ECG data of one pilot in the first group were analyzed over a continuous 35-minute period, during which a total of 2,723 R-peaks were detected. A representative 100-second segment is shown in Fig. 3, clearly demonstrating the accuracy and stability of the detection approach.

ECG index extraction

Eleven heart rate variability (HRV) features were extracted from the ECG signals, encompassing time-domain, frequency-domain, and nonlinear indicators (Table 5). After preprocessing to remove artifacts and ectopic beats, RR intervals were detected using the Pan–Tompkins algorithm²⁹. To enable spectral analysis, the RR interval series were interpolated at a frequency of 4 Hz, which is commonly used and provides adequate resolution to capture both low- and high-frequency components^31,32. Subsequently, frequency-domain features were computed using the Fast Fourier Transform (FFT), while time-domain features were directly derived from the RR series.

In addition, nonlinear features were obtained through Poincaré plot analysis, which quantifies the short- and long-term variability of heart rate by calculating SD1, SD2, and their ratio SD2/SD1. All feature extraction procedures were implemented in MATLAB R2021a using built-in functions and custom scripts to ensure analytical consistency and reproducibility.

Scale data processing

Before the experiment, a survey was conducted to assess the comfort and workload levels of all pilots while in a resting state. All subjects were in a load-free stage during this assessment. The actual scores were calculated by taking the weighted average of the six factors in the scale. The intra-group average scores were then computed to obtain the statistical results, with standard deviation used as the error bars to represent the variability within the data, as shown in Fig. 4.

According to the statistical results of the NASA-TLX scale in two groups of experiments, the workload of flying in foggy days was higher than that of flying in sunny days, with a score of 3.99 ± 0.59 in sunny days and 4.92 ± 0.80 in foggy days, indicating statistical differences between the two groups (t = − 7.198, p < .01). By corresponding the scale score to the flight period, this study discovered that the workload of the pilot increased in the second half of takeoff and cruise and in the landing stage in foggy days. After the experimental videos were reviewed, the high-load physiological characteristics and unsafe driving behaviors of the subjects were analyzed. The result showed that 3/4 of such phenomena occurred in the foggy group, mostly in the third stage of cruising and the landing stage, which was consistent with the statistical data of the scale. Hence, the scale reflected the workload changes of pilots, and the statistical data were reliable.

A 10-point rating system was adopted for each dimension of the NASA-TLX scale to facilitate rapid judgment by pilots during the experiment. The overall workload score was calculated as the arithmetic mean of the six subscales. For example, Subject 01’s pre-takeoff ratings were: mental demand = 3, physical demand = 2, temporal demand = 2, performance = 3, effort = 1, and frustration = 1, resulting in a mean score of 2.17.

To enable segmentation and statistical analysis, workload scores were categorized into three levels. Although the NASA-TLX does not prescribe standard thresholds, we derived empirical cutoffs based on the score distribution in our sample and informed by workload classification logic in prior studies^33,34. Specifically, scores in the range of [0–4] were defined as Level 1 (low workload), (4–7] as Level 2 (moderate workload), and (7–10] as Level 3 (high workload). These thresholds ensured non-overlapping boundaries and practical interpretability.

Based on this scheme, physiological data from 40 pilots were segmented into 60-second samples (N = 30,000; sampling rate = 500 Hz), yielding 453 samples labeled as low workload, 1,135 as moderate, and 512 as high workload. To evaluate differences in physiological responses across workload levels, we applied the Kruskal–Wallis H test³⁵, a nonparametric method suitable for comparing three or more independent groups when data do not meet normality assumptions. When significant effects were found, Dunn’s post-hoc test with Bonferroni correction^36,37 was used to conduct pairwise comparisons and control the family-wise error rate. These analyses supported the selection of discriminative features for subsequent workload modeling.

Screening of evaluation indexes

To examine physiological responses to varying workload conditions, ECG-derived HRV features were segmented into 60-second windows and labeled according to NASA-TLX workload levels. Given that each participant contributed data across all three levels, a repeated-measures design was adopted.

The normality of HRV distributions was assessed using the Shapiro–Wilk test, which indicated significant deviations from normality. Accordingly, the Friedman test was applied to evaluate overall differences in HRV indicators across the workload conditions³⁸. For indicators showing significant omnibus effects, Wilcoxon signed-rank tests were subsequently conducted for pairwise comparisons³⁹. Bonferroni correction was employed to adjust for multiple comparisons and control the family-wise error rate⁴⁰.

To account for inter-individual variability, all HRV features were baseline-corrected using each participant’s resting ECG values collected during a 5-minute pre-task period. All statistical analyses were conducted using IBM SPSS Statistics 26.0.

As a result, four HRV indicators—pNN20, LF/HF, SD2/SD1, and HR—demonstrated significant differences across workload levels. The detailed pairwise test results are reported in Table 6.

Workload analysis

After characteristic indexes were extracted from the ECG data collected in flight simulation, the extracted intragroup ECG indexes per minute of 40 pilots were averaged to exclude individual differences. The experimental comparison results are exhibited in Fig. 5.

Comparison of the overall trends of the two groups of experiment in Fig. 5 showed that the standard deviation of RR interval (SDNN) and HR were affected more significantly by foggy weather than by sunny weather, and the values under the rest state were greater than those in the flight stage. The two values declined and then rose after takeoff. SDNN is an important index reflecting the balance adjustment of the cardiac autonomic nervous system. An increase in SDNN mainly implies an enhanced activity of the sympathetic nervous system and a tension state of the human body. The mean value of HR showed an overall rising trend, and the workload gradually grew, indicating that the pilots’ tension was intensified by the foggy weather, which also increased their workload. This condition was more obvious in the second half of the landing stage and cruising stage and worthy of attention in actual flight.

A significance analysis was performed on the abovementioned ECG indexes for the stages with low and high workloads. Differences were found in LF/HF (t = − 2.246, p < .05), HR (t = − 27.625, p < .01), SD2/SD1 (t = − 7.955, p < .01), and pNN20 (t = − 4.928, p < .01), which could be used to evaluate pilots’ workload.