A method for analyzing time-varying control behavior of a driver based on hilbert-huang transform

Abstract

The application discloses a kind of based on hilbert-huang transform's analysis method of driver time-varying control behavior.Firstly, through improved empirical mode decomposition EMD and causal decomposition method based on EEMD, time domain signal is decomposed into a finite number of IMF, main causal mechanism in driver behavior is effectively extracted.Then, hilbert-huang transform is used to analyze the time-frequency characteristics of signal, and the control behavior of driver is understood in depth through causal analysis.Finally, using strong adaptive model identification method, based on the idea of two-step estimation method, the time-varying driver model parameters are identified, and the identification results are modeled precision evaluation.This method can not only accurately analyze the control behavior of driver, but also effectively process time-varying non-stationary signal, compared with wavelet, fourier transform and other methods, can reach high precision in time and frequency two dimensions simultaneously, help to improve flight safety.

Description

Technical Field

[0001] This invention belongs to the field of driver behavior analysis, and more specifically relates to a method for analyzing time-varying driver control behavior based on Hilbert-Huang transform. Background Technology

[0002] Traditional aircraft pilot behavior analysis primarily relies on fixed configurations or preset patterns (such as cruise, takeoff, and landing). However, this often overlooks a crucial scenario: when an aircraft experiences a sudden change, such as a malfunction, pilots may alter their operational behaviors to adapt to the new flight environment. In such cases, pilot behavior exhibits significant time-varying characteristics. Traditional analytical methods often fail to handle this time-varying behavior because they are mostly based on time-invariant linear or quasi-linear pilot models.

[0003] Furthermore, traditional analysis methods often rely on time-frequency analysis techniques such as Fourier transform or wavelet transform. However, these methods tend to have lower accuracy when processing non-stationary signals and cannot fully extract their frequency characteristics. Therefore, these methods are not ideal for analyzing the time-varying behavior exhibited by drivers.

[0004] In view of the above problems, this invention proposes a novel method for analyzing pilot time-varying control behavior based on the Hilbert-Huang transform. The Hilbert-Huang transform is a very powerful time-frequency analysis tool, particularly effective in handling nonlinear and unstable signal analysis problems. It is an adaptive method that can adapt to changes in signal characteristics at different times. Therefore, using this method to analyze the pilot's time-varying behavior allows for a better understanding of the pilot's operational strategies in the face of aircraft malfunctions, thereby improving flight safety and efficiency. Summary of the Invention

[0005] This invention, through pilot time-varying control behavior analysis and modeling based on Hilbert-Huang transform, can provide a basis for evaluating the flight quality of human-machine systems under time-varying conditions and better design aircraft control laws.

[0006] To achieve the above objectives, the present invention employs the following technical solution: treating the steady-state phase before and after a flight malfunction as a time-invariant system, and the post-malfunction transitional adaptation phase as a time-varying system, the Hilbert-Huang transform method is used to analyze pilot control behavior. The analysis method includes:

[0007] An improvement is made to the Empirical Mode Decomposition (EMD) in the traditional Hilbert-Huang transform.

[0008] The causal decomposition method based on EEMD is adopted. The noise-assisted data analysis method EEMD is used to decompose the time domain signal into a finite number of IMFs and identify the causal interaction of instantaneous phase correlation encoding between two signals at a specific time scale, so as to further analyze the mechanism of driver control behavior.

[0009] Time-frequency characteristic analysis based on Hilbert-Huang, including time-frequency characteristic analysis and causal relationship analysis;

[0010] The time-varying driver control model is identified, the frequency response of the driver's time-varying model is estimated, the McRuer driver model is fitted, the constant parameters of the steady-state segment and the time-varying parameters of the transition adaptation segment are identified in segments, and the modeling accuracy of the identification results is evaluated.

[0011] In one embodiment, the improvements include:

[0012] An improved empirical mode decomposition method, CEEMDAN, is a complete set of empirical mode decompositions with adaptive noise. The key idea of ​​CEEMDAN is to introduce different noises in each round of EMD operation and reduce mode aliasing for each IMF component through multiple iterations. The introduction of adaptive noise helps to better handle uncertainties and noise components in the signal.

[0013] In one embodiment, the causal decomposition method for EEMD includes:

[0014] (1) Add auxiliary white noise to the original bivariate time series signal, and then perform EEMD to decompose a pair of time series A and B into two sets of IMFs and determine the instantaneous phase coherence between each pair of IMFs.

[0015] (2) Remove the IMFs from the given time series A, recombine the process to generate a new set of IMFs (IMF A′), and recalculate the instantaneous phase coherence between the original IMFs (IMF B) and the recombine IMF A′;

[0016] (3) Determine the absolute causal strength and the relative causal strength by calculating the relative ratio of the variance-weighted Euclidean distance between the phase coherence in the original time series and the recombined time series.

[0017] In one approach, the time-frequency characteristic analysis includes:

[0018] First, the time-frequency signal is preprocessed to reduce noise, which facilitates subsequent characteristic analysis and model identification; that is, the CEEMDAN method is used to decompose the acquired raw signal into multiple intrinsic mode functions (IMFs).

[0019] Then, correlation coefficients were calculated for the steady-state segments before and after the fault to obtain the main frequency component of the driver's operation. By analyzing the differences between the main frequency component and its instantaneous energy and the characteristics of the steady-state segments before and after the fault, the driver's transitional adaptation segment to the fault was determined.

[0020] In one approach, the causal analysis utilizes an EEMD-based causal decomposition method to analyze the causal relationship between IMF components of two signals at different time periods and frequencies, thereby gaining a deeper understanding of the driver's control behavior.

[0021] In one approach, the time-varying driver control model identification employs a highly adaptive Hilbert-Huang transform analysis method, based on a two-step estimation approach, to identify the parameters of the time-varying driver model; the specific steps are as follows:

[0022] (1) Based on the driver-in-the-loop simulation results, the Hilbert-Huang transform method is used to estimate the driver frequency response function at different times, i.e., the driver nonparametric model;

[0023] (2) Based on the frequency response recognition results and prior knowledge, select a suitable driver model form; determine the optimization index, select a suitable optimization algorithm to match the time-varying model parameters, and obtain the driver parameter model;

[0024] (3) Apply the driver parameter model for simulation calculation, verify the model by comparing the simulation output with the actual output of the driver, and establish an accuracy evaluation formula to evaluate the identification effect.

[0025] In one approach, the driver's frequency response function uses empirical mode decomposition to adaptively sieve the signal into a series of intrinsic mode function (IMF) components from high to low frequency, selects the component with a higher correlation coefficient as the dominant mode component, and performs a Hilbert transform to obtain the analytical signal.

[0026] In one approach, the objective of optimizing the driver parameter model is to make the amplitude frequency, phase frequency, and frequency response of the driver model identical in the frequency domain, using the minimum weighted sum of squared deviations as the criterion:

[0027]

[0028] Where Φ l , Φ m These represent the frequency response function calculated based on experimental data and the frequency characteristics of the driver model, respectively; A represents the driver's amplitude-frequency characteristic, and P represents the driver's phase-frequency characteristic; (pi / 180) 2 ω is the weighting coefficient between amplitude and phase angle, which aims to ensure that the influence of amplitude and phase angle on the objective function is consistent; ω represents frequency, and the subscript i represents the corresponding frequency point, which is determined by the frequency of the dominant mode component obtained when identifying the frequency response.

[0029] In one approach, establishing an accuracy evaluation formula to assess the identification performance includes:

[0030] The evaluation metrics for time-domain identification accuracy are as follows:

[0031]

[0032] Among them, u human In a human-in-the-loop simulation experiment, u represents the actual control output signal of the driver. model To identify the output signal of the obtained parameter model, var(u human ) represents u human The variance of VAF; when the value of VAF is 1, it is perfect identification. The closer the value is to 1, the smaller the difference between the two signals and the higher the identification accuracy.

[0033] Beneficial effects of this invention:

[0034] This invention presents an analysis method for time-varying driver control behavior based on Hilbert-Huang transform. This method can effectively extract features of non-stationary signals and achieve high analytical accuracy in both time and frequency dimensions, thereby enabling more accurate analysis of driver control behavior.

[0035] By simulating sudden failures in experiments, this invention can effectively analyze and interpret the time-varying control behavior of pilots under sudden aircraft failures and establish a time-varying pilot model, which is of great significance for understanding pilot behavior patterns and providing predictions of future pilot behavior.

[0036] This invention can accurately evaluate the identification effect of a time-varying driver model by identifying its frequency domain response, thereby ensuring the accuracy of the identification model.

[0037] Time-varying pilot control models can provide more accurate assessments of human-machine system flight quality, which helps in designing better aircraft control laws and improving flight safety. Attached Figure Description

[0038] Figure 1 This is a structural diagram of the human-machine closed-loop system for pitch tracking tasks in this invention;

[0039] Figure 2 Set up the schematic diagram for the experimental task;

[0040] Figure 3 This is a schematic diagram showing the locations of the three fault occurrence points;

[0041] Figure 4 For pitch angle tracking curves;

[0042] Figure 5 To track the error signal e;

[0043] Figure 6 The driver input signal u;

[0044] Figure 7 The result of the empirical pattern decomposition;

[0045] Figure 8 For Hilbert's time spectrum;

[0046] Figure 9 This is a three-dimensional time-frequency energy diagram of HHT;

[0047] Figure 10 Here is a flowchart of the CEEMDAN algorithm;

[0048] Figure 11 The modal decomposition result of the control signal u for t*=40s;

[0049] Figure 12 The correlation coefficient for the steady-state period before and after the fault (t* = 40s);

[0050] Figure 13 The instantaneous energy and Fourier spectrum of the main frequency component IMF8 at t*=40s;

[0051] Figure 14 The instantaneous energy and Fourier spectrum of the main frequency component IMF9 at t*=40s;

[0052] Figure 15 The Hilbert time spectrum of the control signal u for t* = 40s;

[0053] Figure 16 This is the result of the separability test;

[0054] Figure 17 The results are from the orthogonality test.

[0055] Figure 18 The causal relationship between e and u before the fault;

[0056] Figure 19 The causal relationship between transition segment e and u;

[0057] Figure 20 To accommodate the causal relationship between e and u;

[0058] Figure 21 The Fourier spectrum of the IMF component obtained from the pre-fault e-decomposition.

[0059] Figure 22 To adapt to the Fourier frequency of the IMF component obtained after e-decomposition;

[0060] Figure 23 Flowchart for driver experiment modeling;

[0061] Figure 24 Flowchart for driver experiment modeling;

[0062] Figure 25 The driver's frequency response every 2 seconds during the experimental transition period;

[0063] Figure 26 The driver gain change is t* = 40s;

[0064] Figure 27 The change in the advance compensation time constant is t*=40s;

[0065] Figure 28 For comparison of the open-loop simulation outputs of the experimental and model drivers;

[0066] Figure 29 This is to compare the closed-loop simulation output of the experimental and model drivers. Detailed Implementation

[0067] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention are described clearly and completely below. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0068] The overall concept of this invention is as follows: Based on the Hilbert-Huang transform time-frequency analysis method, this invention analyzes the adaptive time-varying control behavior characteristics of pilots. Taking a human-in-the-loop experiment under sudden aircraft failure as an example, it analyzes and establishes a time-varying pilot model to verify the effectiveness of the method. Existing studies on pilot control behavior analysis are mostly linear or quasi-linear, belonging to the time-invariant type. However, this method, targeting the study of time-varying behavior, adopts the Hilbert-Huang method, which, compared with wavelet transform and Fourier transform methods, can achieve high accuracy in both time and frequency dimensions, and can more accurately extract the features of non-stationary signals.

[0069] The specific process is as follows: First, a sudden failure is simulated using aircraft configuration switching to design and complete a pilot-in-the-loop simulation experiment. The Hilbert-Huang transform is used for experimental data processing and analysis to obtain the Hilbert time spectrum and the relative causal strength between signals, thus studying the time-varying control behavior characteristics of the pilot. Then, a two-step estimation method based on the Hilbert-Huang transform is used to identify the parameters of the time-varying pilot model, obtaining the pilot's time-varying control model. The pilot model is then placed in open-loop and closed-loop human-machine systems for simulation evaluation to verify the rationality of the pilot model and the effectiveness of the Hilbert-Huang transform method in analyzing the pilot's time-varying control behavior. Through the analysis and modeling of the pilot's time-varying control behavior based on the Hilbert-Huang transform, a basis can be provided for evaluating the flight quality of the human-machine system under time-varying conditions, leading to better design of aircraft control laws.

[0070] Driver-in-the-loop simulation experiment:

[0071] When an aircraft experiences a sudden malfunction, the pilot adaptively changes their control strategy, exhibiting time-varying control behavior. Therefore, a flight simulator is used to conduct a pilot-in-the-loop simulation experiment. By simulating a sudden aircraft malfunction, data on the pilot's time-varying control behavior is collected to support subsequent analysis.

[0072] The experiment was designed as a pitch angle compensation tracking task, with the pilot and aircraft forming a closed-loop human-machine system, such as... Figure 1 As shown, c, e, u, and y represent the pitch angle command input signal, pitch angle tracking error signal, pilot control signal, and aircraft pitch angle output signal, respectively. The flight simulator's head-up display (HUD) shows the error signal between the pitch angle command signal and the aircraft's actual output pitch angle. The pilot makes corresponding maneuvers based on the error signal to eliminate the error. When the aircraft's dynamic characteristics change abruptly, the pilot will also change their control strategy to minimize the error; the pilot's maneuvers exhibit adaptive time-varying characteristics.

[0073] The input signal uses a set of superimposed cosine signals as a multiharmonic function, such as Figure 3 As shown. The analytical expression of this signal is shown in equation (1). From the driver's perspective, it is a random signal that is not easily predicted.

[0074]

[0075] The amplitude A is selected based on the consistency between the power distribution of the multiharmonic signal and the random signal. k and orthogonal frequency ω k =K2π / T, where T is the duration of the experiment (T = 96 s), and the spectral density of the random signal is... Variance σ i 2 =4deg 2 These 15 frequencies ωk It covers a frequency range from 0.26 rad / s to 15.71 rad / s. The selected frequencies and amplitudes are widely used in research at the Moscow Aviation Institute. Its detailed composition is shown in Table 1:

[0076] Table 1. Detailed list of command signal components

[0077] <![CDATA[ω i (rad / s)]]> 0.2618 0.5236 0.7854 1.0472 1.3090 1.5708 2.0944 2.6180 k 9 10 11 12 13 14 15 <![CDATA[ω i (rad / s)]]> 3.1416 3.9270 5.2360 6.2832 7.8540 10.472 15.708

[0078] In the experiment, sudden aircraft failures were simulated by switching aircraft configurations at specified times. The time-varying human-machine system structure is as follows: Figure 2 As shown, when the aircraft configuration changes abruptly, the pilot's control behavior also changes accordingly, demonstrating its nonlinear time-varying characteristics.

[0079] The flight configuration was taken from the Neal-Smith and HAVE PIO configuration libraries. An experiment was designed to simulate a fault in the aircraft control system. This is because the switched configuration adds a first-order inertial element to the flight control system transfer function, resulting in a phase delay that simulates sudden autopilot or flight control computer failures during actual flight. The experimental setup is shown in Table 2.

[0080] Table 2 Fault Simulation Experiment Setup

[0081]

[0082] Set the typical fault occurrence time: t* = 40s, corresponding to the command to increase the pitch angle, such as... Figure 3 As shown.

[0083] Time-domain results analysis:

[0084] The time-domain results of the experiment are as follows Figures 4 to 6 As shown, it includes pitch angle tracking curves, tracking error signals, and driver control signals.

[0085] As can be seen from the figure, both the error signal and the pilot's control signal changed significantly before and after the fault occurred. On the one hand, the sudden change caused by the fault was unpredictable for the pilot. When the sudden change occurred, the pilot had a refractory period and it was difficult to change the control mode in time. At this time, the control input was large, so the tracking error increased sharply at the point of the fault. On the other hand, the change of the aircraft configuration from HP21 to HP25 added a first-order inertial filter element 1 / (s+1), which brought phase lag and reduced the sensitivity of the controlled object. Therefore, it was difficult to effectively track high-frequency signals after the fault, the tracking delay increased, the pilot's control amplitude increased accordingly, and the error also increased.

[0086] Time-frequency analysis of driver behavior based on Hilbert-Huang transform:

[0087] Treating the steady-state phases before and after the fault as a time-invariant system, the driver's control behavior is analyzed using the Hilbert-Huang transform method.

[0088] In control tasks such as aircraft malfunctions and flight mode switching, pilot maneuvers and human-machine interface characteristics exhibit significant time-varying features. Traditional signal analysis methods are mostly based on the stationarity assumption, providing only statistical averages in the time or frequency domains, but failing to simultaneously reveal local features in both. Signals in time-varying systems are often non-stationary; therefore, it is necessary to develop signal analysis methods suitable for extracting time-varying features from non-stationary signals, namely, time-frequency analysis methods. Using this method, the signal is presented in a three-dimensional space of time-frequency-amplitude / energy, revealing both its constituent frequency components and its time-varying characteristics.

[0089] The Hilbert-Huang Transform (HHT) is a relatively new time-frequency analysis method, particularly suitable for analyzing non-stationary and nonlinear signals. Its key feature is the use of Empirical Mode Decomposition (EMD) to effectively decompose complex non-stationary signals into multiple Intrinsic Mode Functions (IMFs), thereby describing the local characteristics of the signal in the time-frequency plane. Compared to other time-frequency analysis methods such as Short-Time Fourier Transform and Wavelet Transform, the HHT method is more adaptive, allowing for the setting of characteristic time scales based on the signal's inherent characteristics. It also avoids the limitations of the Fourier Transform approach and is not bound by the Heisenberg uncertainty principle, thus achieving high accuracy in both time and frequency dimensions and more accurately extracting the features of non-stationary signals.

[0090] The HHT method mainly consists of two parts: Empirical Mode Decomposition (EMD) and Hilbert Transform (HT), with EMD being its core component. Simply put, the basic process of HHT for processing non-stationary signals is as follows: First, the original signal is decomposed into multiple IMF components from high to low frequency using EMD. Then, a Hilbert transform is performed on each IMF component to solve for the instantaneous frequency and instantaneous amplitude, and the corresponding Hilbert time-frequency energy spectrum is plotted.

[0091] 1) Traditional Hilbert-Huang Transform

[0092] The traditional HHT method consists of the following two parts:

[0093] ① Empirical Mode Decomposition (EMD)

[0094] Empirical Mode Decomposition (EMD) is essentially a signal stabilization process. Based on the local mean characteristics and time scale of the original signal, it divides the signal into a series of intrinsic mode function (IMF) components from high to low frequency, such as... Figure 7 As shown.

[0095]

[0096] ② Hilbert Transform (HT)

[0097] HT is a commonly used signal analysis method, essentially defined as the convolution of a signal x(t) with 1 / t, and can emphasize the local properties of x(t), as shown below:

[0098]

[0099] Obtain the analytical signal Where the instantaneous amplitude and instantaneous phase are respectively

[0100]

[0101] Thus, the instantaneous frequency is obtained as

[0102] The instantaneous amplitude and frequency of each IMF component are obtained through HT. The amplitude is then represented on the time-frequency plane to obtain the Hilbert spectrum, which is a three-dimensional graph of the signal's time, frequency, and energy (amplitude). Figure 8 and Figure 9 As shown.

[0103] Essentially, this method performs signal stabilization, aiming to decompose fluctuations or trends with different time scales in the signal step by step to obtain a series of intrinsic mode components. Therefore, the result obtained by performing the Hilbert transform based on these components has certain physical meaning, namely, the distribution law of signal energy on time or spatial scales.

[0104] In summary, analytical methods such as Fourier transform and wavelet transform rely on prior basis functions. In contrast, the HHT method is an adaptive time-frequency localization analysis method that does not require fixed prior basis functions but decomposes the signal based on its inherent characteristics. This makes signal analysis more flexible and versatile, and more suitable for handling non-stationary and nonlinear signals. Furthermore, the intrinsic mode functions can be considered as inherent vibrational modes, and the instantaneous frequencies obtained through the Hilbert transform have clear physical meaning, reflecting certain local characteristics of the signal to some extent. Moreover, the instantaneous frequency is defined as the derivative of the phase function, without considering the entire waveform. Therefore, it is possible to distinguish singular signals from low-frequency signals, which is a significant improvement over methods such as wavelet transform.

[0105] 2) Improved EMD methods

[0106] However, traditional EMD methods may cause mode aliasing problems, that is, the intrinsic mode functions (IMFs) obtained by decomposition contain excessively large characteristic time scales, resulting in multiple oscillation modes in the selected IMFs and multiple frequency components in a single IMF component, thus failing to decompose the signal correctly.

[0107] Therefore, an improved empirical mode decomposition method, Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN), is introduced. This method effectively solves the mode aliasing problem inherent in EMD, improves the stability and accuracy of mode decomposition, and makes it more suitable for the analysis of complex nonlinear and non-stationary signals.

[0108] The key idea of ​​CEEMDAN is to introduce different noise in each round of EMD operation and reduce mode aliasing for each IMF component through multiple iterations. The introduction of adaptive noise helps to better handle uncertainties and noise components in the signal. A schematic diagram of the CEEMDAN algorithm is shown below. Figure 13 As shown.

[0109] The steps of the CEEMDAN algorithm are as follows:

[0110] Let E i Let be the i-th IMF component obtained after EMD decomposition, and let be the i-th IMF component obtained after CEEMDAN decomposition. z(t) is the original signal to be decomposed, v j Let be Gaussian white noise with a mean of 0, j = 1, 2, ..., N, where is the number of times the white noise is added, and ε is the weighting coefficient of the Gaussian white noise. Adding paired positive and negative Gaussian white noise to the original signal z(t) yields the new signal z(t)+(-1). q εv i (t), where q = 1 or 2. Let the first-order eigenmode function obtained by EMD of the new signal be C1:

[0111]

[0112] (1) The average of the generated N intrinsic mode components is obtained by summing:

[0113]

[0114] (2) Calculate the residual after removing the first mean mode component:

[0115]

[0116] (3) Gaussian white noise is then added to the residual r1(t) to obtain another new signal. EMD decomposition is then performed on this signal to obtain...

[0117] First-order eigenmode function From this, we can obtain the second-order modal components:

[0118]

[0119] (4) Calculate the residual after removing the second mean mode component:

[0120]

[0121] (5) Repeat the above steps until the obtained residual signal is a monotonic function and cannot be further decomposed, at which point the algorithm ends.

[0122] If the number of intrinsic mode components obtained at this point is K, then the original signal z(t) is decomposed into:

[0123]

[0124] Compared to EEMD and CEEMD, CEEMDAN introduces noise adaptively and reduces mode aliasing through iterative EMD, which helps reduce noise residue and effectively improves the purity of the decomposition results. The improvements to CEEMDAN make it more suitable for processing complex nonlinear and non-stationary signals. For signals that traditional methods may struggle to handle, CEEMDAN provides a more effective decomposition method.

[0125] 3) Causal decomposition method based on EEMD

[0126] Existing methods for detecting causal relationships in time series are primarily based on Bayesian prediction. However, causes and effects are often simultaneous, and most real-world causal interactions are reciprocal, such as predator-prey relationships and physiological regulation of bodily functions.

[0127] The causal decomposition is based on two assumptions: (1) any causal relationship between the source signal and the target signal can be quantified by instantaneous phase dependence and decomposed into intrinsic components at a specific time scale; (2) the phase dynamics in the target signal generated by the source signal are separable from the target itself.

[0128] Among them, EEMD is a noise-assisted data analysis method used to further improve the separability of IMFs during the decomposition process, and defines the true IMF components as the average value of the test set.

[0129]

[0130] In the above formula, the additional noise level r determines the separability of the IMF components. Therefore, it is necessary to determine the appropriate noise level coefficient r through separability and orthogonality tests. In principle, the separability should be maximized (the root mean square of the pairwise correlation between IMFs should be minimized (<0.05)) while maintaining acceptable non-orthogonal leakage (<0.05).

[0131] Based on Galileo's covariant principle of causality, the causal relationship between two time series is defined as follows: the cause is the input, and the effect is the output; removing the cause results in no effect. Therefore, if the intrinsic component of B that is causally related to A is removed from B itself, the instantaneous phase correlation between A and B decreases, then variable A causes variable B, and vice versa. This decomposition and recombination process can quantify the differential causal relationship between the corresponding IMFs of two time series.

[0132]

[0133] Coh(A,B′) <Coh(A,B)~Coh(A′,B) (13)

[0134] The specific steps of this method are as follows:

[0135] (1) Add auxiliary white noise to the original bivariate time series signal, and then perform EEMD to decompose a pair of time series A and B into two sets of IMFs and determine the instantaneous phase coherence between each pair of IMFs.

[0136] (2) Remove the IMFs from the given time series A, recombine the process to generate a new set of IMFs (IMF A′), and recalculate the instantaneous phase coherence between the original IMFs (IMF B) and the recombine IMF A′;

[0137] (3) Determine the absolute causality intensity D(S) by calculating the relative ratio of the variance-weighted Euclidean distance between the phase coherence in the original time series (IMFA, IMF B) and the recombined time series (e.g., IMFA′, IMF B). 1j →S 2j ), D(S 2j →S 1j ) and relative causal strength C(S) 1j →S 2j ), C(S 2j →S 1j ).

[0138]

[0139] Therefore, this application introduces this analytical method, using EEMD to decompose the time-domain signal into a finite number of IMFs, and identifies the causal interaction of the instantaneous phase correlation encoding between two signals at a specific time scale, further analyzing the mechanism of driver control behavior. Based on this, by decomposing the error signal e and the driver control signal u, the differences in the causal correlation between the IMF components of the two signals in the time-varying and time-invariant phases can be analyzed, as well as the driver's tracking effect on different frequency components of the signal in different phases.

[0140] Time-frequency characteristic analysis based on Hilbert-Huang:

[0141] Next, the experimental data will be further processed and analyzed based on the Hilbert-Huang time-frequency method:

[0142] (1) Time-frequency characteristic analysis

[0143] First, in order to improve the quality and accuracy of the signal, the above signal needs to be preprocessed for noise reduction, which will facilitate subsequent characteristic analysis and model identification.

[0144] The so-called noise reduction preprocessing is to use the CEEMDAN method to decompose the acquired raw signal into multiple intrinsic mode functions (IMFs). Since the higher the noise content, the lower the correlation coefficient, the correlation coefficient can be used to remove the modal components with excessive noise, thereby achieving the noise reduction effect.

[0145] The correlation coefficients between each modal component and the original signal are calculated using the following formula:

[0146]

[0147] In the formula, x i Represents the intrinsic modal components, y i This represents the original signal u or e. Since the correlation coefficient is applicable to the calculation of stationary signals, the original signal needs to be segmented. Assuming the transition adaptation period is 10 seconds after the fault occurs, the steady-state periods before and after the fault are calculated separately.

[0148] Taking the control signal u obtained with a switching time of t*=40s as an example, the Empirical Mode Decomposition (EMD) of the signal is performed using the CEEMDAN method, and a series of IMF components are obtained by sieving the signal from high to low frequency. Figure 11 As shown.

[0149] Then, correlation coefficients were calculated for the steady-state periods before and after the fault, as shown below. Figure 12 The results are shown below. Among them, the IMF components with the highest correlation coefficients, namely IMF8 and IMF9, are considered as the dominant frequency components of the driver's maneuver. Their time-domain instantaneous energy change curves and Fourier spectra are plotted, as shown below. Figure 13 and Figure 14As shown, the main frequency components before and after the fault are 2.91 rad / s and 1.57 rad / s, respectively.

[0150] On the one hand, the instantaneous energy of the higher frequency component IMF8 (approximately 2.91 rad / s) drops sharply at the point of failure, indicating that the configuration switch significantly impacts the pilot's control. The large-delay configuration HP25 makes it extremely difficult for the pilot to track the higher frequency components. On the other hand, the energy of the lower frequency component IMF9 (1.57 rad / s) surges at the point of failure, and this frequency band further increases in the subsequent steady-state phase, suggesting that the pilot focuses more on tracking the low-frequency components in the latter half of the test.

[0151] Taking both into account, such as Figure 13 and Figure 14 As shown, the IMF component and its instantaneous energy within the red box (approximately 40-50 seconds) are different from the characteristics of the preceding and following steady-state segments. Therefore, these 10 seconds can be regarded as the driver's transitional adaptation segment to the fault in this case, and the above assumption is reasonable.

[0152] Furthermore, the correlation coefficients of the IMF8-12 components are all greater than 0.2, and they can be considered as dominant mode components for subsequent Hilbert time-frequency spectrum plotting (see...). Figure 15 Other components were discarded due to their low correlation coefficients and high noise levels. As can be seen from the figure, after a fault occurs, high-frequency energy weakens while low-frequency energy strengthens. This means that due to response delay, the driver has difficulty effectively tracking high-frequency signals and therefore focuses more on tracking the low-frequency band.

[0153] (2) Causal relationship analysis

[0154] Using the EEMD-based causal decomposition method, we can analyze the causal relationship between IMF components of two signals at different time periods and frequencies, which can be understood as correlation, thereby gaining a deeper understanding of the driver's control behavior. The following example, taking the case of t* = 40s in Experiment 1, demonstrates a causal relationship analysis between the driver's input signal (error e) and the output signal u.

[0155] Because the e and u signals contain a lot of noise, to ensure the separability of the IMF components, it is necessary to determine the appropriate noise level coefficient r through separability and orthogonality tests. The principle is to maximize separability (minimize the root mean square of the pairwise correlation between IMFs (<0.05)) while maintaining acceptable non-orthogonal leakage (<0.05). Taking the steady-state period before the fault (1-40s) as an example, we obtain... Figure 16 and Figure 17 Based on the test results, the noise level coefficient r should be selected as 0.65.

[0156] Causal decomposition between e and u is performed in three stages: before the fault (1-40s), during the transition (40-50s), and after adaptation (50-96s). The first eight IMF components, not lower than the fundamental frequency of the command signal, are extracted to obtain the corresponding relative casual strength (RCS). Figures 18 to 20 As shown, blue bars represent the relative causal strength of u to e, and red bars represent the opposite. Specific values ​​for the relative causal strength of u to e are shown in Table 3. A ratio of 0.5 indicates no causal relationship was detected, or that there was no difference in causal strength when the two signals were mutually causal. A ratio close to 0 or 1 indicates a strong causal influence from signal e or signal u, respectively.

[0157] In addition, the Fourier spectra corresponding to the IMF5-IMF8 components in each steady-state segment are calculated as follows: Figure 21 and Figure 22 As shown in the figure. Combining the relative causal strength histogram and Table 3, it can be seen that, on the one hand, before the fault occurred, the driver had the best tracking effect on the frequency component IMF5 of about 3.14 rad / s, with an RCS as high as 0.93; while after adaptation, the driver had the best tracking effect on the frequency components IMF6 and IMF7 of 1.53 rad / s and 0.99 rad / s, respectively, with a decrease in the dominant frequency, which is consistent with the fact that the configuration delay increased after the fault, making it difficult for the driver to track high frequencies.

[0158] On the other hand, the causal strength of the transition adaptation phase is significantly lower than that of the steady-state phase, indicating that the nonlinearity and non-stationarity of this phase are relatively high, making it difficult for the driver to maintain good tracking, especially the higher frequency components of the signal.

[0159] Table 3. Relative causal strength of u to e

[0160] 1-40s 0.93 0.60 0.65 0.26 40-50s 0.50 0.51 0.20 0.76 50-96s 0.51 0.77 0.23 0.56

[0161] Identification (modeling) of time-varying driver control model

[0162] Based on the analysis results of Hilbert-Huang's time-varying control behavior of drivers, the frequency response of the time-varying driver model is estimated, the McRuer driver model is fitted, the constant parameters of the steady-state segment and the time-varying parameters of the transitional adaptation segment are identified in segments, and the modeling accuracy of the identification results is evaluated.

[0163] Hilbert-Huang-based time-varying driver model identification method:

[0164] This application employs a highly adaptive Hilbert-Huang transform analysis method, based on a two-step estimation approach, to identify the parameters of a time-varying driver model. The specific steps are as follows:

[0165] 1. Based on the driver-in-the-loop simulation results, the Hilbert-Huang transform method is used to estimate the driver's frequency response function at different times, i.e., the driver's nonparametric model;

[0166] 2. Based on the frequency response recognition results and prior knowledge, select a suitable driver model form; determine the optimization index, select a suitable optimization algorithm to fit the time-varying model parameters, and obtain the driver parameter model;

[0167] 3. Simulation calculations are performed using a driver parameter model. The model is validated by comparing the simulation output with the driver's actual output. An accuracy evaluation formula is then established to assess the recognition effect.

[0168] Driver experimental modeling process as follows Figure 23 As shown.

[0169] (1) Identify frequency response

[0170] Since the power of the input function is limited to specific input frequencies, the driver's response function can only be identified at these frequencies. Therefore, combining the principle of the Hilbert-Huang transform, empirical mode decomposition is used to adaptively sieve the signal into a series of intrinsic mode function (IMF) components from high to low frequency. The component with the higher correlation coefficient is selected as the dominant mode component, and a Hilbert transform is applied to it to obtain the analytic signal. Analogous to wavelet time-frequency analysis, under the ideal condition that the input function is a multi-harmonic function, the time-varying driver frequency response function can be approximated as the ratio of the Hilbert-Huang transforms of the control signal u and the error signal e:

[0171]

[0172] Since the analytic signal is a complex number, and the transfer function is also a complex number, it can be further written in the following form:

[0173]

[0174] Where a and These correspond to the amplitude and phase of the analytical signal, respectively. Based on this, a time-varying nonparametric model of the driver can be obtained, that is, the driver's frequency response at different times.

[0175] (2) Model Form Selection

[0176] Human drivers are nonlinear biological systems. However, in continuous control tasks, when properly trained and given constant environmental conditions, the driver's manual control behavior can be represented by a quasi-linear response function H. P (This relates the perceived tracking error to the control input u) and is described by the residual signal n (all control inputs unrelated to the command signal) representing nonlinear behavior, i.e., the McRuer quasi-linear driver model (see...). Figure 24 This model has become a powerful research tool for studying the impact of different perceived information on pilot control behavior, evaluating aircraft handling quality, and assessing the design of different control systems. Furthermore, the McRuer quasi-linear pilot model has fewer time-varying parameters, which facilitates the study of time-varying pilot control behavior. Therefore, this application will apply this model for parameter identification.

[0177] Figure 24 The response function H in P The format is as follows:

[0178] H p (s,t)=K p (t)(1+T L (t)s)e -τs H nm (19)

[0179]

[0180] Among them, K P T represents driver gain. L The driver's advance compensation time constant represents both of these parameters, and the neuromuscular system H... nm The driver's reaction time delay τ takes into account the limitations of the driver's manipulation and physiological conditions. To reduce the number of unknown parameters in the driver model and facilitate the subsequent identification of time-varying parameters, the neuromuscular system and delay parameters are set to constant values, and ω is selected. nm =10 rad / s, τ = 0.2s, and the model ignores the lag element. Further verification will be conducted to confirm that the accuracy of the model still meets the requirements.

[0181] (3) Identify model parameters

[0182] After determining the form of the driver model, the driver model parameters for different times can be fitted based on the estimated frequency response. Before optimizing the model parameters, appropriate optimization indices need to be determined. The goal of driver parameter model optimization is to make the amplitude frequency, phase frequency, and frequency response of the driver model close to the frequency response in the frequency domain. Intuitively, the index that minimizes the weighted sum of squared deviations can be used as follows:

[0183]

[0184] Where Φ l , Φ m These represent the frequency response function calculated based on experimental data and the frequency characteristics of the driver model, respectively; A represents the driver's amplitude-frequency characteristic, and P represents the driver's phase-frequency characteristic; (pi / 180) 2ω is the weighting coefficient between amplitude and phase angle, which aims to ensure that the influence of amplitude and phase angle on the objective function is consistent; ω represents frequency, and the subscript i represents the corresponding frequency point, which is determined by the frequency of the dominant mode component obtained when identifying the frequency response.

[0185] (4) Parameter identification accuracy index

[0186] After identifying the driver model parameters, it is necessary to define an identification accuracy index to evaluate the parameter identification effect. On one hand, open-loop simulation is performed using the same error input signal, comparing the control signal obtained from the time-varying driver parameter model with the driver's actual control output. On the other hand, the time-varying driver parameter model is placed in a closed-loop human-machine system for simulation, and the simulation results are compared with human-machine simulation experimental data. The time-domain identification accuracy evaluation index is as follows:

[0187]

[0188] Among them, u human In a human-in-the-loop simulation experiment, u represents the actual control output signal of the driver. model To identify the output signal of the obtained parameter model, var(u human ) represents u human The variance of VAF. A VAF value of 1 indicates perfect identification; the closer the value is to 1, the smaller the difference between the two signals, and the higher the identification accuracy. This index characterizes the accuracy with which the identified parameter model describes the driver's actual control behavior. For quasi-linear models, it includes both the accuracy of parameter fitting and the linearity of the driver's control behavior.

[0189] Time-varying driver model identification results:

[0190] Based on the two-step estimation method, the driver's frequency response at different times is first identified to obtain the time-varying nonparametric model results. Then, the results are used to fit the McRuer quasi-linear model to identify the time-varying driver model parameters under sudden faults.

[0191] (1) Nonparametric model identification results:

[0192] The series of IMF components obtained by EMD correspond to different frequencies. The IMF components with high correlation to the original signal are selected and subjected to Hilbert transformation to obtain the amplitude and phase that change with time. Then, the driver frequency response that changes with time is obtained according to Equation (17), that is, the nonparametric model is identified.

[0193] like Figure 25 The figure shows the driver's frequency response curves every 2 seconds during the transition phase. The amplitude-frequency response changes little at different times, but the phase-frequency response changes significantly. The driver's frequency response characteristics will be represented in the driver's mathematical model by the time-varying parameter K.P and T L It is reflected in.

[0194] (2) Parameter model identification results

[0195] Using the obtained frequency response to fit the mathematical model of the driver, the time-varying driver model parameters under sudden fault conditions are identified, namely the driver gain K. P and advance compensation time constant T L Since the driver's continuous operation can be divided into a steady-state segment and a transitional adaptation segment, the steady-state signal is considered as a statistically quasi-stationary signal, and K... P and T L The parameters were treated as constants for fitting, while during the transition adaptation period, they were treated as time-varying parameters and fitted at each time step. For the group with t* = 40s, the transition adaptation period was set to 10s.

[0196] On the one hand, constant fitting was performed on the two steady-state segments before and after the fault, and the model parameter identification results are shown in Table 4. On the other hand, for the driver's transition adaptation segment, time-varying parameters were identified with a step size of 0.1s, and the moving average method was used to analyze parameter K. P and T L The identification results are smoothed to obtain the final model parameter change curves for the entire task phase.

[0197] Table 4. Results of steady-state model parameter identification in the experiment.

[0198]

[0199] A comparison of the steady-state periods before and after the fault is shown in Table 4. P and T L Both showed significant improvements, indicating that the delay caused by the control system fault was substantial, necessitating a simultaneous increase in K. P and T L To improve the speed of system response.

[0200] During the transitional adaptation phase, such as Figure 26 and Figure 27 As shown, the driver gain K P The trend shows an initial increase followed by a decrease. In the early stages of the malfunction, aircraft delay increases, requiring pilots to increase K... P To accelerate the system's response speed, K needs to be appropriately reduced in the later stages of the transition phase. P To reduce system overshoot. And T L The trend is that it first increases and then decreases, with T increasing in the early stages of the fault. L This can reduce overshoot, overcome oscillations, and accelerate the system's response speed; as the system gradually stabilizes in the later stages of the transition, T can be appropriately reduced. L This increases the ability to suppress disturbances.

[0201] Therefore, the time-varying control behavior of the driver under abrupt faults can be approximated as a proportional-derivative controller with time-varying parameters, where the driver gain K is... P and advance compensation time constant T L This reflects the driver's control behavior characteristics. For system delays caused by malfunctions, the driver increases K... P and T L To speed up the system's response time, but K P An excessively large value for T may result in significant overshoot, leading to oscillations. L Excessive stress can cause disturbances and make the system more sensitive, and the control difficulty varies at different fault locations. Therefore, during the transitional adaptation phase, the driver needs to continuously adjust K. P and T L The size is adjusted to speed up response and maintain system stability.

[0202] (3) Validation of time-varying driver parameter model

[0203] On the one hand, since the input to the open-loop system is the experimentally obtained error signal, which is the same as the input used when identifying the model, it is not affected by feedback loops. Therefore, open-loop simulation can better reflect the identification accuracy of the time-varying transfer function model at each moment. On the other hand, human-machine closed-loop system simulation can verify the overall identification effect of the driver model. Therefore, this application places the driver model in both human-machine open-loop and human-machine closed-loop systems for simulation evaluation.

[0204] Open-loop and closed-loop simulation results are as follows Figures 28 to 29 As shown in the figure, the trends of driver input changes in the experimental and identified models are basically consistent, indicating good identification performance. However, the identification performance for the flat portion of the actual driver control signal is poor. This is because the driver's control strategy assumes that there is no need to actively track the error signal when it fluctuates around 0°, and this part belongs to the nonlinear residual in the driver's control behavior.

[0205] The VAF values ​​for the open-loop and closed-loop simulation identification accuracy of the above experiments are shown in Table 5. It can be seen that, due to the presence of closed-loop feedback, the accuracy of the closed-loop simulation is slightly lower than that of the open-loop simulation, reaching above 0.65 in all stages. In contrast, the accuracy of the open-loop simulation reaches above 0.7 in all stages. This indicates that the identification results of the time-varying parameters of the driver model are relatively accurate, which also demonstrates the effectiveness of the Hilbert-Huang transform analysis of the driver's time-varying control behavior.

[0206] Table 5. Model parameter identification accuracy index (VAF value) for the experiment.

[0207] Open loop 0.7859 0.8654 0.8560 closed loop 0.7093 0.8225 0.8083

[0208] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for analyzing driver time-varying control behavior based on Hilbert-Huang transform, characterized in that: Treating the steady-state period before and after the flight malfunction as a time-invariant system, the pilot's control behavior is analyzed using the Hilbert-Huang transform method. The analysis method includes: An improvement to the traditional Empirical Mode Decomposition (EMD) in the Hilbert-Huang transform; The causal decomposition method based on EEMD is adopted. The noise-assisted data analysis method EEMD is used to decompose the time domain signal into a finite number of IMFs and identify the causal interaction of instantaneous phase correlation encoding between two signals at a specific time scale, so as to further analyze the mechanism of driver control behavior. Time-frequency characteristic analysis and causal relationship analysis based on Hilbert-Huang transform; The time-varying driver control model is identified, the frequency response of the driver's time-varying model is estimated, the McRuer driver model is fitted, the constant parameters of the steady-state segment and the time-varying parameters of the transition adaptation segment are identified in segments, and the modeling accuracy of the identification results is evaluated. The time-varying driver control model identification employs a highly adaptive Hilbert-Huang transform analysis method, based on a two-step estimation approach, to identify the parameters of the time-varying driver model; the specific steps are as follows: (1) Based on the driver-in-the-loop simulation results, the Hilbert-Huang transform method is used to estimate the driver frequency response function at different times, i.e., the driver nonparametric model; (2) Based on the frequency response recognition results and prior knowledge, select a suitable driver model form; determine the optimization index, select a suitable optimization algorithm to match the time-varying model parameters, and obtain the driver parameter model; (3) Apply the driver parameter model for simulation calculation, verify the model by comparing the simulation output with the actual output of the driver, and establish an accuracy evaluation formula to evaluate the identification effect; The driver frequency response function uses empirical mode decomposition to adaptively sieve the signal into a series of intrinsic mode function (IMF) components from high to low frequency, selects the component with higher correlation coefficient as the dominant mode component, and performs Hilbert transform to obtain the analytical signal. The goal of optimizing the driver parameter model is to make the amplitude frequency, phase frequency, and frequency response of the driver model identical in the frequency domain. The criterion used is to minimize the weighted sum of squared deviations, as follows: ； in , These are the frequency response function and the frequency characteristics of the driver model, respectively, calculated based on experimental data. Represents the driver's amplitude-frequency characteristics. This represents the driver's phase frequency characteristics; The weighting coefficient between amplitude and phase angle is used to ensure that the influence of amplitude and phase angle on the objective function is consistent. Represents frequency, subscript The corresponding frequency point is determined by the frequency of the dominant mode component obtained when identifying the frequency response. The establishment of an accuracy evaluation formula to assess the recognition effect includes: The evaluation metrics for time-domain identification accuracy are as follows: ； in, The actual control output signal of the driver in a human-in-the-loop simulation experiment. To identify the output signal of the obtained parameter model, represent The variance of VAF; when the value of VAF is 1, it is perfect identification. The closer the value is to 1, the smaller the difference between the two signals and the higher the identification accuracy.

2. The method for analyzing time-varying driver control behavior based on Hilbert-Huang transform according to claim 1, characterized in that: The improvements include: An improved empirical mode decomposition method, CEEMDAN, is a complete set of empirical mode decompositions with adaptive noise. The key idea of ​​CEEMDAN is to introduce different noises in each round of EMD operation and reduce mode aliasing for each IMF component through multiple iterations. The introduction of adaptive noise helps to better handle uncertainties and noise components in the signal.

3. The method for analyzing time-varying driver control behavior based on Hilbert-Huang transform according to claim 1, characterized in that: The causal decomposition method of EEMD includes: (1) Add auxiliary white noise to the original bivariate time series signal, and then perform EEMD to decompose a pair of time series A and B into two sets of IMFs and determine the instantaneous phase coherence between each pair of IMFs; (2) Remove the IMFs from the given time series A, recombine them to generate a new set of IMFs A′, and recalculate the instantaneous phase coherence between the original IMF B and the recombine IMF A′; (3) Determine the absolute causal strength and the relative causal strength by calculating the relative ratio of the variance-weighted Euclidean distance between the phase coherence in the original time series and the recombined time series.

4. The method for analyzing time-varying driver control behavior based on Hilbert-Huang transform according to claim 1, characterized in that: The aforementioned time-frequency characteristic analysis includes: First, the time-frequency signal is preprocessed to reduce noise, which facilitates subsequent characteristic analysis and model identification; that is, the CEEMDAN method is used to decompose the acquired raw signal into multiple intrinsic mode functions (IMFs). Then, the correlation coefficients were calculated for the steady-state periods before and after the fault.

5. The method for analyzing time-varying driver control behavior based on Hilbert-Huang transform according to claim 1, characterized in that: The aforementioned causal analysis utilizes the EEMD-based causal decomposition method to analyze the causal relationship between IMF components of two signals at different time periods and frequencies, thereby gaining a deeper understanding of the driver's control behavior.

Images