Research Article  Open Access
On Flight Risk Quantitative Evaluation under Icing Conditions
Abstract
The quantitative assessment of flight risk under icing conditions was taken as the research object. Based on multifactor coupling modeling idea, the pilotaircraftenvironment coupling system was built. Considering the physical characteristics and randomness of aircraft icing, the extreme values of critical flight risk parameters were extracted by the Monte Carlo flight simulation experiment. The flight characteristics were studied comprehensively and heavytail characteristics and the distributions of different flight parameters were verified. Flight risk criterion was developed and onedimensional extreme flight risk probability was calculated. Further, in order to solve the limitation of onedimensional extreme value, with the Copula theory, the joint distribution model of flight parameters with three distinct distribution types was built. The optimal Copula model was selected by identification of unknown parameters and goodness of fit tests, and the threedimensional extreme flight risk probability was defined. Based on the quantitative flight risk, the accident induction mechanism under icing conditions was discussed. Airspeed and roll angle under asymmetry icing conditions were more sensitive and had a more significant impact on flight safety. This method can provide reference for safety manipulation, boundary protection, and risk warning during icing flight.
1. Introduction
Despite the significant improvement in flight safety for the advanced flybywire aircraft, the risk of encountering upset conditions remains a serious threaten. Aircraft icing is a critical external environmental factor, which can induce a lossofcontrol incident if the pilot does not respond in the correct manner [1]. The United States National Transport Safety Board (NTSB) had aircraft icing on its “Most Wanted” list of safety recommendations from 1997 [2]. In the most serious icing accident on June 1, 2009, Air France’s A330 airliner crashed over the Atlantic Ocean, and all 228 people on board were killed. The cause of the accident was ice formation that caused the autopilot to turn off, and the pilot’s mishandling led to a stall. On February 11, 2018, a Russian An148 plane crashed shortly after takeoff and all 71 people on board got killed, due to aircraft icing. How to objectively evaluate flight risk under random icing conditions is a fundamental and difficult issue.
In terms of icing contamination on flight safety, the current researches focus mainly on the following.
(1) Icing Simulation Method. By the icing wind tunnel experiment [3–5], numerical simulation [6–8], and flight test [9, 10], the icing similarity criterion, icing meteorological conditions and formation mechanism, and ice shape evolution process [11–14] were studied.
(2) Effect of Icing on Aerodynamic and Flight Mechanics Characteristics. The change rules of aerodynamic parameters and stability derivatives and nonlinear and stall conditions after icing were studied [15–17]. The dynamic response characteristics under different ice severity and distribution conditions [18, 19] and icing influence on the overall performance, maneuverability and stability of the aircraft [20] were explored.
(3) Icing Tolerance Flight and Boundary Protection. The parameters identification method with aircraft icing [21, 22], the ice detection and warning method [23], the manipulation method and control law design after icing [24, 25], and the ice management system [26, 27] were hot research issues.
(4) AntiIcing System and Ice Airspace Avoidance. The efficient and effective antiicing methods [28], accurate prediction, and avoidance meteorology methods of icing airspace [29] were intensively studied and applied to engineering practice. Related institutions increased the rules and procedures. For example, some relevant regulations were added to the icing airworthiness review in FAA 25.
The above research and literature have high academic value, but there are still few relevant theories and methods in the evaluation and prediction of flight risk probability under icing conditions. The dynamic assessment of flight risk under icing conditions is at the centre of current aircraft safety and airworthiness work, and it is an urgent problem to improve the existing flight safety analysis theory. The existing safety regulations and guidelines, like SAE ARP4761 [30], SAE ARP4754A [31], MILHDBK516B [32], and MILSTD882E [33], have clear ideas and methods for assessing the probability of flight accidents due to the failure of the internal hardware system, but lack quantitative indicators for the dynamic flight risk induced by the external environment. The current research in the theory and method of flight risk assessment under the influence of the external environment is mostly qualitative analysis or static reliability assessment. Such methods have obvious limitations in assessing the flight risk probability under icing conditions: (1) based on the deterministic model, they cannot reflect the randomness and uncertainty of the multifactor coupling system under icing conditions; (2) these methods cannot evaluate the dynamic flight risk based on the flight process and difficult to estimate and predict the quantitative probability index for the flight risk caused by icing; (3) flying with ice is a small probability, high risk flight subject, which is impossible to get enough samples data under the same actual conditions. The sample size required for the quantitative analysis of the risk probability is about 3.81 × 10^{7}3.85 × 10^{8} (taking 95% confidence levels). How to assess the dynamic flight risk objectively and quantitatively based on limited samples is an urgent problem to be resolved.
On the other hand, Copula theory has become a research hotspot for multivariate extreme problems and has been applied in small probability events’ evaluation and prediction such as finance, disaster early warning, and signal processing [34–36]. Because Copula function can construct flexible multivariate distribution function and describe the correlation among different marginal distribution parameters, especially nonlinear, asymmetric, and tail correlation, the multivariate extreme Copula model is much suitable for assessing events with high hazards and small probability.
Based on the distributed simulation method [37] under icing conditions and the flight risk evaluation method [8] of two flight parameters with the same type mathematical distributions, we further proposed a quantitative flight risk assessment method of several flight parameters with different mathematical distributions in the case of multifactor coupling. The proposed method could intuitively indicate the critical factors restricting flight safety and flight risk grade. In view of these, the organization of the paper is as follows. In Section 2, in order to obtain the flight state parameters whose variation could directly reflect the movement state and flight risk, the pilotaircraftenvironment realtime simulation platform is constructed. Then, the extremum of flight parameters is calculated by the Monte Carlo simulation and flight safety critical parameters are defined. And also, the credibility of the simulation method is verified by comparing the simulation results undertaken by the real pilot in the loop with the pilot model produced by the Monte Carlo simulations under the same conditions. The flight safety critical parameters’ distribution characteristics are analyzed and the flight risks based on onedimensional parameters are calculated in Section 4. In Section 5, in order to emerge from the limitations of onedimensional extreme value model, the construction method of multivariate extreme Copula model under icing conditions is explored and the flight risk probability quantitatively is calculated. Finally, the conclusion is drawn in Section 6.
2. Modeling of Closed Loop Simulation System under Icing Conditions
Flying with ice involves the coupling characteristics of the longitudinal and lateral and other nonlinear characteristics. To accurately simulate the dynamic response of the aircraft, the sixdegreeoffreedom equations are built based on the F16 model structure [38] and the wind tunnel experimental data of a certain transporter model. The actuator model can be simplified as a firstorder lag with rate and position limits, which were validated to respond correctly and in accordance with the Transport Class Model (TCM) [39]. And also, control augmentation is adopted to improve the control characteristics. Close response characteristic between the two models demonstrates the effectiveness of the established aircraft simulation system.
2.1. Aircraft Dynamic Model
The nonlinear dynamic model of the aircraft can be expressed as
where is state vector with airspeed, angle of attack, sideslip angle, quaternion, roll rate, pitch rate, yaw rate, and space position parameters; namely, . is the control vector, including throttle deflection, elevator deflection, aileron deflection, and rudder deflection; namely, .
In order to avoid singularities in the calculation process, the quaternion method is used to construct the aircraft dynamic model, as shown in the following.
where are the components of the velocity in the body coordinate system, respectively; are, respectively, the component of the resultant force in the body coordinate system of aircraft; are the component of the resultant moment in the body coordinate system; () are the intermediate variables related to the moment of inertia.
The attitude angle of the aircraft: can be obtained by quaternion transformation and the sum of squares of quaternions is equal to 1.
2.2. Pilot Model
Pilot manipulation is essential for flight safety; therefore, it is of importance for the pilot to perform correctly and timely. The transfer function pilot model proposed by McRuer has been widely used as inwhere , , , , and are pilot behavior characteristic parameters. represents the pilot’s static gain, which is about 1100; represents the pilot’s delayed response time to sense the change of the aircraft attitude, which is about 0.060.3 s and related to nerve conduction and stimulation. For multidegree manipulation, the value is larger than 0.2 s. represents the pilot’s lead compensation time constant, which is about 01s and reflects the pilot’s prediction of the dynamic response; represents the time delay constant for the transmission and processing of the central information, reflecting the load on the pilot, about 01 s; indicates muscle nerve lag time, about 0.1 ± 20% s. The different values for the parameters in (7) reflect the randomness of different pilot manipulations. How are the values to represent the randomness identified; in other words, what is the distribution of the pilot manipulation parameters? Statistics of the actual pilot operation behavior show that follows the lognormal distribution and the other four parameters obey the truncated normal distribution. Figures 1 and 2 show the distribution of and .
Under icing conditions, the pilot manipulation could be divided into three stages: the pilot senses the change of flight state and begins to manipulate after seconds; pilot handles the stick, throttle, and pedal to modify the altitude until the kinematic parameters begin to reverse, here time comes to and then keep the attitude of the aircraft stable. Pilot input and output models under icing conditions can be expressed as (8) and Figure 3, where and are different McRuer models based on the database of pilot compensation manipulation behavior.
2.3. Icing Effect Model
Through lots of icing experimental researches on DHC6, Bragg [26] developed a physically representative model to analyze the flight dynamics of an iced aircraft, as shown in the following.where is an arbitrary performance, stability, or control derivative affected by ice accretion; represents the severity of icing, and the greater , the greater the effect of ice on the aerodynamic parameters. Usually the value of is about 00.3. is a constant for a certain aircraft, which is usually obtained by flight test or flight simulation calculations.
In actual flight, considering the randomness of ice formation and shedding, there will be asymmetrical life and drag produced by twoside wings, and further additional roll moments and yaw moments are generated, which is particularly problematic and dangerous when oneside deicing system fails. Lampton [20] developed an asymmetric icing effect model to simulate the difference in force and moment, as shown in where and are lift coefficient and drag coefficient affected by aircraft icing; is the distance along the body Xaxis from the mean geometric chord to the aircraft centerline; is the dynamic pressure and is the wing area.
Then (10) can be integrated into the 6DOF equations to simulate the flight dynamic characteristic of the aircraft under asymmetrical icing conditions. What needs to be emphasized is that aircraft icing not only changes the aerodynamic parameters, but also deteriorates allowable range of some critical safetyrelated parameters. Here, take stall angle of attack (AOA) as an example and the model [20] is shown as Here, is stall angle of attack affected by aircraft icing.
It is necessary to note that the above icing effect models are simplified to some extent and suitable for preliminary analysis of the ice accumulation effect on aerodynamic parameters. If highprecision numerical simulation is necessary, the flight safety parameter ranges under different ice severities and flight states should be recorded through icing wind tunnel tests or flight tests, stored in the form of a database, and called by interpolation. The data in this paper were obtained by CFD simulation and compared with the wind tunnel test data, which showed the simulation was credible. As this is not the focus of this article, it will not be discussed in detail.
3. Monte Carlo Simulation Process and Credibility Analysis
3.1. Monte Carlo Simulation Process
The current general method of predicting flight risk is tantamount to observe the variation of some critical safety parameters, because flight risk is usually accompanied by abnormal flight parameters. Considering the randomness of pilot model, external environment models, and icing effect model, combined with the nonlinear simulation system, the flight of a typical icing condition is simulated and flight parameter extreme values are extracted based on the Monte Carlo method. The simulation process is illustrated in Figure 4.
(1) Simulation Status Setting Module. Set the initial flight state of the aircraft, the total number of simulation (), icing time, and the typical distribution of aircraft icing.
(2) Random Variable Extraction Module. The external environment and pilot maneuvering parameters required for the th simulation are extracted based on the Monte Carlo method. According to a certain area’s external environment database, the aircraft ice formation and external disturbance factors (such as a sudden wind, turbulence) could be determined, and the corresponding aerodynamic parameters could be extracted based on the icing database. Of course, if we need to study the dynamic effects of icing environment on flight simulation more accurately, we need to consider the impact of icing dynamic process on flight aerodynamics, such as the dynamic crosslinking of icing CFD simulation and aerodynamic derivatives.
(3) Flight Simulation Module. Based on the pilot model, icing effect model, environment model, and the aircraft model, the flight simulation platform is established with the RungeKutta algorithm. The simulation step length is 0.02s. The flight dynamics model and flight control system are loaded on the realtime simulation computer constructed by LinksTech and connect with other simulation knots by the reflective memory network.
(4) Extreme Parameters Storage Module. In each flight simulation, the maximum or minimum value of the status parameters is stored in the extreme value database, such as the maximum AOA and the minimum airspeed during the th simulation.
(5) Simulation Process Control Module. Let , and determine whether is greater than the setting number . If not, the simulation returns to the second step. Experiment shows that when is greater than 2000, the statistical characteristic is stable. So, set equal to 2500 to compromise operation time and result stability.
3.2. Safety Critical Parameters Definition
Select an initial state of level flight with flight height of 3000 m, airspeed 120 m/s. Aircraft enters icing airspace after 5 s and the right deicing system is failed. The variations of flight parameters in the simulation are shown in Figure 5.
(a)
(b)
(c)
(d)
(e)
(f)
From Figure 5, with the ice accumulation on the wings and before the pilot intervenes, the velocity drops slightly due to the lift decreases and drag increases; the aircraft rolls to the right side and begins to yaw, due to the asymmetrical lift and drag caused by the failure of the right deicing equipment; and at the same time, flight height rapid decreases. Then, the pilot senses the change of flight state and intervenes to manipulate. The thrust increases first to obtain energy and then decreases to try to recover the initial state. During the process, the AOA increased by 103.3%, the rolling angle reached 49.9 degrees, the speed margin reduced by 58.3%, and the height loss by 116.3 meters. During the same time, the stall AOA decreased, the lateral control hinge torque increased, and the pilot was more difficult to navigate. Here, we choose the AOA, roll angle, and airspeed as the flight extremum parameters, which drastic change and have a significant impact on flight safety. Of course, if the pilot performs a takeoff or landing mission, the height loss is particularly critical. For different missions, different flight parameters are selected to represent the flight risk.
3.3. Monte Carlo Simulation Credibility Analysis
About the limits of natural icing conditions and the limits of statistical stability for the flight extreme parameters, it is impossible to obtain extreme values through more than 2000 flight tests under the same natural conditions. Also, it is difficult to perform more than 2000 tests using the pilotintheloop ground simulation platform. It is feasible to obtain the parameter extremum through the Monte Carlo simulation method. Judging the reliability of this method is the primary problem to be resolved.
200 tests are performed, respectively, with the real pilot in the loop and the pilot model produced by the Monte Carlo simulations under the same conditions to verify the credibility [40]. The distribution of the extracted extreme parameters for the two different kinds of experiments is compared in Figure 6, which shows that the distribution is similar. To further and more clearly verify the similarity distributions, the QuantilesQuantiles (QQ) plots, Rsquare method, correlation coefficient method, and KolmogorovSmirnov (KS) test are performed as Figure 7 and Table 1 show.

The QQ graphical method plots the quantile scatter of the twosample empirical distribution function on the same graph. The straight lines in Figure 7 indicate that the two samples belong to the similar distribution. KS test is to determine the goodness of fit of the two samples by calculating the empirical probability difference. The P value of the KS test represents the test statistic obtained from the observation of the original sample and the minimum significance level at which the original hypothesis can be rejected. All the P values are greater than 0.05, which means there is generally no reason to reject the original hypothesis. The coefficient of correlation and Rsquare between the samples is greater than 0.9, indicating that the twosample data showed a strong linear correlation. So it is possible to pass the test and accept the original hypothesis. Considering comprehensively, it can be concluded that the pilot model simulation data has the same distribution type as the real test data of the pilot in the loop, and the pilot model simulation data sample can be used to evaluate the flight risk probability.
4. Distribution Characteristics of OneDimensional Flight Extremum Parameters
4.1. Statistical Characteristics Analysis of Flight Extremum Parameters
The Monte Carlo simulation method was used to extract the extreme value samples of flight parameters. At a certain confidence level, the statistical characteristics of the extremum parameters were analyzed, as shown in Table 2, in order to study the distribution characteristics of the flight extreme parameters.

In Table 2, the kurtosis coefficients above 3 show that the extremum samples of the three parameters are more concentrated and have longer tails than the normal distribution; the skewness coefficients show that the left tail of the minimum velocity extremum is longer and the right tails of the roll angle extremum and the AOA extremum are longer. The three flight parameters have heavy tails and different distribution types.
4.2. Mathematical Distribution Type Analysis
This obvious heavytailed distribution is common in lowfrequency and high risk events. For this type of heavytailed distribution, the most effective description method is extreme value theory [35, 36], which can model the tail of the extreme value probability distribution of a random sequence, and calculate the specific overlimit probability value based on the mathematical model. In extreme value theory, extreme value (EV) distribution, general extreme value (GEV) distribution can describe the heavytailed distribution. Lognormal (Logn) distribution, Weibull distribution, and exponential (Exp) distribution can also describe the heavytailed characteristics. Here, five typical distribution families are selected to identify onedimensional extreme value parameters and the normal distribution is also listed to compare, as shown in the following.The identified weight parameters for the above equations are listed in Table 3.

In order to judge the accuracy of the above models, the goodness of fit test is applied by the KS test, AndersonDarling (AD) test, and χ^{2} test. The results are shown in Tables 4–6.



In Table 4, KS test values of EV and Weibull distributions are both less than 0.1, and P values are both greater than 0.05, indicating that 95% confidence levels can pass the test, and distribution fitting degree is higher. Furthermore, AD test is more sensitive to heavytailed distribution, Weibull and EV distribution have smaller test values, and Weibull’s P value is greater than 0.05. Therefore, by integrating various test rules, the identification accuracy of Weibull distribution is higher. Similarly, the accuracy of GEV distribution is the highest in Tables 5 and 6. In summary, the minimum velocity extremum accords with the Weibull distribution, and the rolling angle and the AOA extremum meet the GEV distribution.
4.3. Flight Risk Probability Based on OneDimensional Extreme Parameter
Flight risk events are often accompanied by flight parameters exceeding the limit. According to the mathematics distribution of the above three parameters, a criterion for the occurrence of flight risk events is proposed. For the range of these three parameters is quite different, the flight risk criterion is constructed by using the method of normalization of extremum parameters.
By consulting the flight manual of the certain aircraft, the critical value of the AOA is a function of the Mach and the position of the flap. Therefore, according to the aerodynamic data in the icing and failure conditions, and the Mach and flap position when extracting the AOA extreme value, the critical value of AOA is obtained by interpolation calculation, and then the normalized AOA is obtained by . Similarly, the critical value of the roll angle is 85 deg. Take the airspeed corresponding to the maximum lift coefficient in the current flight environment as the critical value. So, the criterion for flight risk under icing condition can be described as
Based on the mathematical models and the identified values of the three flight parameters, the flight risk can be calculated, respectively, by (19)(21). For the above simulation, the probabilities of flight risk are 0.0556, 0.0268, and 0.0051. According to MILSTD882D, the flight risk calculated, respectively, by minimum airspeed and maximum roll angle extreme parameters reaches the B level, namely, “possible”, while, the flight risk calculated by maximum AOA parameter reaches the C level, namely, “occasionally”. The flight risk probabilities based on the different flight parameters in the same flight are different, which means using a single extremum to access flight risk is not comprehensive and it is difficult to choose a suitable parameter. At the same time, we can conclude that, in the asymmetric wing icing conditions, pilot should pay much attention to the airspeed in case of stall and trying to keep lateral balance in case of oversize rolling.
To comprehensive and accurate flight risk evaluation, in the next section, a risk assessment model based on the threedimensional flight extremum parameters with different mathematical distributions is built.
5. ThreeDimensional Extreme Risk Model
The single variable flight risk assessment model has some limitations and can not fully reflect the risk situation. Because the joint distribution of multivariate variables is not only related to the edge distribution of each component, but also to the correlation between variables, the conclusion of the univariate can not be extended directly to multivariate cases. The extreme value character of a single component may not contain the joint extreme value of the entire vector. Copula theory is proposed for the evaluation of extremum distribution and has good applicability for describing the correlation among multiple parameters. This section focuses on the construction of multivariate flight parameter extreme risk model with different mathematical distributions.
5.1. Mathematical Distribution Type Hypothesis
Assume the distribution function of the extremum vector is and the marginal distribution functions are , , . We have concluded that obeys Weibull distribution and and obey GEV distribution at 95% confidence level in Section 4. According to Sklar theorem [8], for any , there must be a Copula function , which satisfies
A threedimensional Copula model is constructed based on asymmetric Archimedes Copula, as shown in where is Archimedes Copula generating function; is marginal distribution of the sample . Firstly, Copula is generated by and ; secondly, Copula is generated by combining and . Thirdly, assuming is the parameter of Copula and is the parameter of Copula , expressions for and are given; finally, according to (23), is brought into . Then, the threedimensional Copula is constructed to calculate flight risk probability. By analogy, more parameters can be considered for comprehensive evaluation. If there are too many parameters, convergence and precocity of Copula weight parameters calculation may be a problem. Furthermore, in the actual, we only need to select the most sensitive parameters for flight safety
Here, , , are all continuous distribution functions, so Copula is the only one. The commonly used multidimensional Copula models and corresponding generators are shown as (24)(28). Gumbel Copula: Frank Copula: Clayton Copula: GS Copula: Joe Copula:
5.2. Identification of the Parameters in Copula Functions
Using nonlinear identification algorithms, the identified weight parameters and for the above models are calculated and listed in Table 7.

Goodness of fit test is applied by the Akaike Information Criteria (AIC), Bayesian Information Criterion (BIC), χ^{2} test, and KS test, in order to judge the Copula model accuracy. The results are shown in Table 8. We can analyze that the test values of Clayton model and Joe model are higher than those of significance levels 0.01, 0.02, and 0.05; that is, the two Copula models could pass the test at 99% confidence level, while other models could not pass the test. Continuing to compare AIC values, BIC values, and chisquare statistics, Joe model can more accurately describe the threedimensional extreme value distribution model under the condition than Clayton model. What is more, Joe model is sensitive to heavytail characteristic. So we choose Joe model to describe the joint probability distribution of minimum airspeed extreme, maximum roll angle extreme, and maximum AOA extreme.

5.3. ThreeDimensional Extreme Parameter Flight Risk Probability and Analysis
According to the risk discriminant (18) and Joe Copula model, the flight risk probability of threedimensional extreme parameters is shown in
represents Joe Copula model based on formula (28). According to the identification results in Table 7, flight risk probability under the above simulation conditions (asymmetric wing icing) is 0.0924. By comparing the flight risk probability calculated by onedimensional extreme GEV model and Weibul model, it can be found that the risk probability calculated by the multivariate extreme Copula model is larger. It shows that the flight risk situation when the extreme values of three different flight parameters exceed the limits is fully considered. Compared with onedimensional extreme value of flight parameters to determine flight risk, multivariate extreme value Copula method is more comprehensive and accurate.
It should be emphasized that the flight risk probability calculated here is similar to the accident rate in SAE ARP4761 [30], which is largely a reference value. Flight risks in different situations can be compared and analyzed horizontally, such as in different harsh environment conditions, in different hardware failure conditions, or in different risk task subjects.
In order to intuitively analyze the relationship between flight risk and extremum parameters, Figure 8 shows the joint probability density distribution of Joe Copula model and Clayton Copula model when =0.8. The maps show strong coupling and high risk probability near the upper tail and low risk probability in other regions. Compared with Clayton Copula model, Joe Copula model has more centralized risk probability distribution and its probability density gradient is larger near the upper tail. The distribution reflects that when multiple flight parameters are faced with overlimit, the flight risk probability is large, while when a single flight parameter is close to overlimit, the flight risk gradient is larger, which requires pilots to precisely and carefully manipulate to avoid cliffbreaking flight risk. From another dimension, Joe Copula model shows that the certain aircraft has a high design level and the flight risk events may occur only when multiple flight parameters exceed the limits with a high probability, which is similar to the flight test results of the certain aircraft.
6. Conclusions
(1) Considering the randomness and uncertainty of the environment and pilot’s manipulation under icing conditions, the feasibility of the Monte Carlo method for flight dynamics simulation under icing conditions is verified and relevant experiments are carried out on a realtime simulation platform.
(2) Extremum samples of key flight risk parameters are extracted through a large number of simulations, and the distribution types of different extremum parameters are obtained. The minimum velocity extremum obeys Weibull distribution, maximum roll angle, and AOA extremum obey GEV distribution. The flight risk based on single extremum parameter is calculated, and the risk values are different. So, the limitation of single extremum parameter in evaluating flight risk is analyzed. Furthermore, combining with the theory of multivariate extremum, the flight risk model based on Joe Copula model is constructed and the correlation between flight extremum parameters is further analyzed.
(3) Through comparative analysis, the minimum speed extremum has a greater impact on the flight risk under icing conditions, and the pilots should pay special attention to the airspeed. When encountering asymmetric icing, pilots need be aware of the abnormal roll of the aircraft, avoid entering upset state, and control the aircraft quickly and smoothly to restore normal attitude.
(4) Flight risk probability values proposed in this paper can be used as an effective supplement to the current flight safety regulations and have positive significance for the review of flight safety and airworthiness. At the same time, because the flight risk probability is a reference value, the flight risk degree under different icing conditions or other failure conditions can be compared and analyzed horizontally. This paper mainly considers the influence of icing on aerodynamic parameters, but does not consider the risk factors of icing on engine, sensor and so on. Therefore, the next step will take into account the flight risk and the warning method under the condition of multiple factors crosslinking after icing, in order to improve the pilot’s situational awareness.
Nomenclature
:  Throttle deflection 
:  Elevator deflection 
:  Aileron deflection 
:  Rudder deflection 
:  Airspeed, m/s 
:  Angle of attack, rad 
:  Sideslip angle, rad 
:  Roll angle, rad 
:  Quaternion, 
:  Roll rate, pitch rate, yaw rate, body axis, rad/s 
:  Space position, m 
:  Components of the velocity, body axis, m/s 
:  Components of the resultant force, body axis, N 
:  Components of the resultant moment, body axis, N·m 
:  The pilot’s static gain 
:  The pilot’s delayed response time, s 
, :  Lift and drag coefficients affected by aircraft icing 
:  Distance along the body Xaxis from the mean geometric chord to the aircraft centerline, m 
:  Dynamic pressure 
:  Wing area, 
:  Flight risk probability. 
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors extend their special thanks to China Flight Test Establishment for the opportunity to test initial ideas. This study was cosupported by the National Natural Science Foundation of China (61503406, 61873351) and the National Key Basic Research Program (2015CB755800).
References
 Y. Cao, W. Tan, and Z. Wu, “Aircraft icing: an ongoing threat to aviation safety,” Aerospace Science and Technology, vol. 75, pp. 353–385, 2018. View at: Publisher Site  Google Scholar
 A. Reehorst, J. Chung, M. Potapczuk, and Y. Choo, “Study of icing effects on performance and controllability of an accident aircraft,” Journal of Aircraft, vol. 37, no. 2, pp. 253–259, 2000. View at: Publisher Site  Google Scholar
 G. E. Fujiwara and M. B. Bragg, “Method for designing hybrid airfoils for icing windtunnel tests,” Journal of Aircraft, vol. 56, no. 1, pp. 137–149, 2019. View at: Publisher Site  Google Scholar
 M. Leng, S. Chang, Y. Lian, and H. Wu, “Experimental study of water film dynamics under wind shear and gravity,” AIAA Journal, vol. 56, no. 5, pp. 1803–1811, 2018. View at: Publisher Site  Google Scholar
 G. E. Fujiwara and M. B. Bragg, “Method for designing threedimensional swept hybrid wings for icing windtunnel tests,” Journal of Aircraft, pp. 1–17, 2018. View at: Publisher Site  Google Scholar
 M. F. Alam, D. S. Thompson, and D. K. Walters, “Hybrid Reynoldsaveraged NavierStokes/largeeddy simulation models for flow around an iced wing,” Journal of Aircraft, vol. 52, no. 1, pp. 244–256, 2015. View at: Publisher Site  Google Scholar
 X. Tong, D. Thompson, Q. Arnoldus, E. Collins, and E. Luke, “Threedimensional surface evolution and mesh deformation for aircraft icing applications,” Journal of Aircraft, vol. 54, no. 3, pp. 1047–1063, 2017. View at: Publisher Site  Google Scholar
 J. M. Wang, H. J. Xu, Y. Xue et al., “Flight risk evaluation of tailplane icing based on extreme value theory,” Acta Aeronaut Astronaut Sincial, vol. 37, no. 10, pp. 3011–3022, 2016. View at: Google Scholar
 T. Ratvasky, B. Barnhart, S. Lee, and J. Cooper, “Flight testing an iced business jet for flight simulation model validation,” in Proceedings of the 45th AIAA Aerospace Sciences Meeting and Exhibit, Aerospace Sciences Meetings, Reno, Nev, USA, 2007. View at: Publisher Site  Google Scholar
 E. Whalen and M. B. Bragg, “Aircraft characterization in icing using flighttest data,” Journal of Aircraft, vol. 42, no. 3, pp. 792–794, 2005. View at: Publisher Site  Google Scholar
 X. Yi, Numerical Computation of Aircraft Icing and Study on Icing Test Scaling Law, China Aerodynamics Research and Development Center Graduate School, 2007.
 T. Hu, H. Lv, B. Tian, and D. Su, “Choosing critical ice shapes on airfoil surface for the icing certification of aircraft,” Procedia Engineering, vol. 80, pp. 456–466, 2014. View at: Publisher Site  Google Scholar
 S. Kumar and E. Loth, “Aerodynamic simulations of airfoils with uppersurface iceshapes,” Journal of Aircraft, vol. 38, no. 2, pp. 285–295, 2001. View at: Publisher Site  Google Scholar
 X. Li, J. Bai, J. Hua, K. Wang, and Y. Zhang, “A spongy icing model for aircraft icing,” Chinese Journal of Aeronautics, vol. 27, no. 1, pp. 40–51, 2014. View at: Publisher Site  Google Scholar
 T. P. Ratvasky, B. P. Barnhart, and S. Lee, “Current methods modeling and simulating icing effects on aircraft performance, stability, control,” Journal of Aircraft, vol. 47, no. 1, pp. 201–211, 2010. View at: Publisher Site  Google Scholar
 L. Qu, Y. Li, H. Xu, D. Zhang, and G. Yuan, “Aircraft nonlinear stability analysis and multidimensional stability region estimation under icing conditions,” Chinese Journal of Aeronautics, vol. 30, no. 3, pp. 976–982, 2017. View at: Publisher Site  Google Scholar
 Y. H. Cao, Z. L. Wu, Y. Su et al., “Aircraft flight characteristics in icing conditions,” Progress in Aerospace Sciences, vol. 74, pp. 62–80, 2015. View at: Google Scholar
 M. Bragg, T. Hutchison, and J. Merret, “Effect of ice accretion on aircraft flight dynamics,” in Proceedings of the 38th Aerospace Sciences Meeting and Exhibit, Aerospace Sciences Meetings, AIAA 20000360, Reno, NV, USA, 2000. View at: Publisher Site  Google Scholar
 H. Kim and M. Bragg, “Effects of leadingedge ice accretion geometry on airfoil performance,” in Proceedings of the 17th Applied Aerodynamics Conference, Fluid Dynamics and Colocated Conferences, AIAA993150, Norfolk, VA, USA, 1999. View at: Publisher Site  Google Scholar
 A. Lampton and J. Valasek, “Prediction of icing effects on the lateral/directional stability and control of light airplanes,” Aerospace Science and Technology, vol. 23, no. 1, pp. 305–311, 2012. View at: Publisher Site  Google Scholar
 R. Aykan, C. Hajiyev, and F. Çalişkan, “Kalman filter and neural networkbased icing identification applied to A340 aircraft dynamics,” Aircraft Engineering and Aerospace Technology, vol. 77, no. 1, pp. 23–33, 2005. View at: Publisher Site  Google Scholar
 R. Aykan, C. Hadjiyev, and F. Caliskan, “Aircraft icing detection, identification and reconfigurable control based on kalman filtering and neural networks,” in Proceedings of the AIAA Atmospheric Flight Mechanics Conference and Exhibit, AIAA20056220, San Francisco, Calif, USA, 2005. View at: Publisher Site  Google Scholar
 F. Caliskan, R. Aykan, and C. Hajiyev, “Aircraft icing detection, identification, and reconfigurable control based on Kalman filtering and neural networks,” Journal of Aerospace Engineering, vol. 21, no. 2, pp. 51–60, 2008. View at: Publisher Site  Google Scholar
 Z. Li, H. Xu, Y. Xue, and B. Pei, “Study on flight safety manipulation space under complex conditions,” Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, vol. 233, no. 2, pp. 725–735, 2019. View at: Publisher Site  Google Scholar
 W. Chen, H. J. Xu, X. L. Wang et al., “Reconfigurable control methods of icing aircraft longitudinal motion based on robust servo LQR,” Acta Aeronauticaet Astrinautica Sinica, vol. 38, no. 1, pp. 71–80, 2017. View at: Google Scholar
 M. B. Bragg, T. Basar, W. R. Perkins et al., “Smart icing systems for aircraft icing safety,” SAE Technical Paper 2003012100, 2003. View at: Google Scholar
 B. Martos, R. Ranaudo, B. Norton et al., “Development, implementation, and pilot evaluation of a modeldriven envelope protection system to mitigate the hazard of inflight ice contamination on a twinengine commuter aircraft,” NASA/CR2014218320, 2014. View at: Google Scholar
 Mengl Wu I, Chiyu Wang, Yunpeng Li, and Q Nie I, “Research on the heating of deicing fluid in a new reshaped coiled tube,” Mathematical Problems in Engineering, vol. 2017, Article ID 3254631, 9 pages, 2017. View at: Publisher Site  Google Scholar  MathSciNet
 B. Pei, H. Xu, and Y. Xue, “Lyapunov based estimation of flight stability boundary under icing conditions,” Mathematical Problems in Engineering, vol. 2017, Article ID 6901894, 10 pages, 2017. View at: Publisher Site  Google Scholar  MathSciNet
 SAE ARP 4761. “Guidelines and methods for conducting the safety assessment process on civil airborne systems and equipment.” Warrendale, PA, USA. 1996.
 SAE ARP 4754A. “Guidelines for development of civil aircraft and systems,” Warrendale, PA, USA. 2010.
 MILHDBK516B, Department of Defense Handbook: Airworthiness Certification Criteria, U.S. Department of Defense, 2005.
 MILSTD882E, Standard practice for system safety, U.S. Department of Defense, 2012.
 Q. Chen, D. Wang, and Mi. Pan, “Multivariate timeVarying GH copula Garch model and its application in the financial market risk measurement,” Mathematical Problems in Engineering, vol. 2015, Article ID 286014, 9 pages, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 T. Berger, “Forecasting valueatrisk using time varying copulas and EVT return distributions,” International Economics, vol. 133, pp. 93–106, 2013. View at: Publisher Site  Google Scholar
 D. Forcellini, “A new methodology to assess indirect losses in bridges subjected to multiple hazards,” Emerging Science Journal, vol. 2, no. 6, p. 400, 2018. View at: Publisher Site  Google Scholar
 Z. Li, H. J. Xu, Y. Xue et al., “Distributed simulation method for manmachineenvironment complex system under icing condition,” Journal of Systems Engineering and Electronics, vol. 40, no. 5, pp. 1167–1174, 2018. View at: Google Scholar
 L. Sonneveldt, Nonlinear F16 model description, University of Technology, Delft, The Netherlands, 2010.
 R. M. Hueschen, “Development of the transport class model (TCM) aircraft simulation from a subscale generic transport model (GTM) simulation,” NASA/TM2011217169, 2011. View at: Google Scholar
 J. B. Osei, M. AdomAsamoah, A. A. Awadallah Ahmed, and E. B. Antwi, “Monte Carlo based seismic hazard model for Southern Ghana,” Civil Engineering Journal, vol. 4, no. 7, p. 1510, 2018. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Zhe Li et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.