# On Flight Risk Quantitative Evaluation under Icing Conditions

Authors: Zhe Li; Haojun Xu; Yuan Xue; Zhe Zhang; Xiaocong Duan
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The quantitative assessment of flight risk under icing conditions was taken as the research object. Based on multifactor coupling modeling idea, the pilot-aircraft-environment 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 heavy-tail characteristics and the distributions of different flight parameters were verified. Flight risk criterion was developed and one-dimensional extreme flight risk probability was calculated. Further, in order to solve the limitation of one-dimensional 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 three-dimensional 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.


---

## Body

## 1. Introduction

Despite the significant improvement in flight safety for the advanced fly-by-wire aircraft, the risk of encountering upset conditions remains a serious threaten. Aircraft icing is a critical external environmental factor, which can induce a loss-of-control 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 An-148 plane crashed shortly after take-off 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) Anti-Icing System and Ice Airspace Avoidance. The efficient and effective anti-icing 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 ARP-4761 [30], SAE ARP-4754A [31], MIL-HDBK-516B [32], and MIL-STD-882E [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 × 107-3.85 × 108 (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 pilot-aircraft-environment real-time 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 one-dimensional parameters are calculated in Section 4. In Section 5, in order to emerge from the limitations of one-dimensional 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 six-degree-of-freedom equations are built based on the F-16 model structure [38] and the wind tunnel experimental data of a certain transporter model. The actuator model can be simplified as a first-order 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(1)x˙=fx,uwherex is state vector with airspeed, angle of attack, sideslip angle, quaternion, roll rate, pitch rate, yaw rate, and space position parameters; namely, x=[V,α,β,q0,q1,q2,q3,p,q,r,xg,yg,zg]T. u is the control vector, including throttle deflection, elevator deflection, aileron deflection, and rudder deflection; namely, u=[δth,δe,δa,δr]T.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.(2)f1=uu˙+vv˙+ww˙Vf2=uw˙-wu˙u2+w2f3=v˙V-vV˙V2cos⁡β(3)f4f5f6f7=-120pqr-p0-rq-qr0-p-r-qp0·q0q1q2q3(4)f8=c1r+c2pq+c3L-+c4N-f9=c5pr-c6p2-r2+c7M-f10=c8p-c2rq+c4L-+c9N-(5)xgygzg=q02+q12-q22-q322q1q2-q0q32q1q3+q0q22q1q2+q0q3q02-q12+q22-q322q2q3-q0q12q1q3-q0q22q2q3+q0q1q02-q12-q22+q32·uvw(6)u˙=rv-qw+Xmv˙=pw-ru+Ymw˙=qu-pv+Zmwhereu,v,w are the components of the velocity in the body coordinate system, respectively; X,Y,Z are, respectively, the component of the resultant force in the body coordinate system of aircraft; L-,M-,N- are the component of the resultant moment in the body coordinate system; ci (i=1,2,…,9) 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 in(7)Gs=Kpe-τsTLs+1TIs+1TNs+1where Kp, τ, TL, TI, and TN are pilot behavior characteristic parameters. Kp represents the pilot’s static gain, which is about 1-100; τ represents the pilot’s delayed response time to sense the change of the aircraft attitude, which is about 0.06-0.3 s and related to nerve conduction and stimulation. For multidegree manipulation, the value is larger than 0.2 s. TL represents the pilot’s lead compensation time constant, which is about 0-1s and reflects the pilot’s prediction of the dynamic response; TI represents the time delay constant for the transmission and processing of the central information, reflecting the load on the pilot, about 0-1 s; TN 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 Kp follows the lognormal distribution and the other four parameters obey the truncated normal distribution. Figures 1 and 2 show the distribution of Kp and τ.Figure 1
Gain (Kp) distribution of the pilot.Figure 2
Time delay (τ) distribution of the pilot.Under icing conditions, the pilot manipulation could be divided into three stages: the pilot senses the change of flight state and begins to manipulate aftert1 seconds; pilot handles the stick, throttle, and pedal to modify the altitude until the kinematic parameters begin to reverse, here time comes to t2 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 F(s) and G(s) are different McRuer models based on the database of pilot compensation manipulation behavior.(8)Us=0tp≤t1Fs·Rs-Yst1<tp≤t2Gs·Rs-Ystp>t2Figure 3
Schematic diagram of pilot manipulation model.


### 2.3. Icing Effect Model

Through lots of icing experimental researches on DHC-6, Bragg [26] developed a physically representative model to analyze the flight dynamics of an iced aircraft, as shown in the following.(9)CAiced=1+ηKCACAwhere C(A) 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 0-0.3. KC(A) 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 two-side wings, and further additional roll moments and yaw moments are generated, which is particularly problematic and dangerous when one-side deicing system fails. Lampton [20] developed an asymmetric icing effect model to simulate the difference in force and moment, as shown in (10)ΔCLice=12CL-12CLiceΔCDice=12CDice-12CDΔLice=QSwΔCLicedmgcΔNice=QSwΔCDicedmgcwhere CLice and CDice are lift coefficient and drag coefficient affected by aircraft icing; dmgc is the distance along the body X-axis from the mean geometric chord to the aircraft centerline; Q is the dynamic pressure and Sw is the wing area.Then (10) can be integrated into the 6-DOF 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 safety-related parameters. Here, take stall angle of attack (AOA) as an example and the model [20] is shown as (11)ΔCL=CLice-CLαstall=fΔCLHere, αstall 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 high-precision 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.


## 2.1. Aircraft Dynamic Model

The nonlinear dynamic model of the aircraft can be expressed as(1)x˙=fx,uwherex is state vector with airspeed, angle of attack, sideslip angle, quaternion, roll rate, pitch rate, yaw rate, and space position parameters; namely, x=[V,α,β,q0,q1,q2,q3,p,q,r,xg,yg,zg]T. u is the control vector, including throttle deflection, elevator deflection, aileron deflection, and rudder deflection; namely, u=[δth,δe,δa,δr]T.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.(2)f1=uu˙+vv˙+ww˙Vf2=uw˙-wu˙u2+w2f3=v˙V-vV˙V2cos⁡β(3)f4f5f6f7=-120pqr-p0-rq-qr0-p-r-qp0·q0q1q2q3(4)f8=c1r+c2pq+c3L-+c4N-f9=c5pr-c6p2-r2+c7M-f10=c8p-c2rq+c4L-+c9N-(5)xgygzg=q02+q12-q22-q322q1q2-q0q32q1q3+q0q22q1q2+q0q3q02-q12+q22-q322q2q3-q0q12q1q3-q0q22q2q3+q0q1q02-q12-q22+q32·uvw(6)u˙=rv-qw+Xmv˙=pw-ru+Ymw˙=qu-pv+Zmwhereu,v,w are the components of the velocity in the body coordinate system, respectively; X,Y,Z are, respectively, the component of the resultant force in the body coordinate system of aircraft; L-,M-,N- are the component of the resultant moment in the body coordinate system; ci (i=1,2,…,9) 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 in(7)Gs=Kpe-τsTLs+1TIs+1TNs+1where Kp, τ, TL, TI, and TN are pilot behavior characteristic parameters. Kp represents the pilot’s static gain, which is about 1-100; τ represents the pilot’s delayed response time to sense the change of the aircraft attitude, which is about 0.06-0.3 s and related to nerve conduction and stimulation. For multidegree manipulation, the value is larger than 0.2 s. TL represents the pilot’s lead compensation time constant, which is about 0-1s and reflects the pilot’s prediction of the dynamic response; TI represents the time delay constant for the transmission and processing of the central information, reflecting the load on the pilot, about 0-1 s; TN 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 Kp follows the lognormal distribution and the other four parameters obey the truncated normal distribution. Figures 1 and 2 show the distribution of Kp and τ.Figure 1
Gain (Kp) distribution of the pilot.Figure 2
Time delay (τ) distribution of the pilot.Under icing conditions, the pilot manipulation could be divided into three stages: the pilot senses the change of flight state and begins to manipulate aftert1 seconds; pilot handles the stick, throttle, and pedal to modify the altitude until the kinematic parameters begin to reverse, here time comes to t2 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 F(s) and G(s) are different McRuer models based on the database of pilot compensation manipulation behavior.(8)Us=0tp≤t1Fs·Rs-Yst1<tp≤t2Gs·Rs-Ystp>t2Figure 3
Schematic diagram of pilot manipulation model.


## 2.3. Icing Effect Model

Through lots of icing experimental researches on DHC-6, Bragg [26] developed a physically representative model to analyze the flight dynamics of an iced aircraft, as shown in the following.(9)CAiced=1+ηKCACAwhere C(A) 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 0-0.3. KC(A) 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 two-side wings, and further additional roll moments and yaw moments are generated, which is particularly problematic and dangerous when one-side deicing system fails. Lampton [20] developed an asymmetric icing effect model to simulate the difference in force and moment, as shown in (10)ΔCLice=12CL-12CLiceΔCDice=12CDice-12CDΔLice=QSwΔCLicedmgcΔNice=QSwΔCDicedmgcwhere CLice and CDice are lift coefficient and drag coefficient affected by aircraft icing; dmgc is the distance along the body X-axis from the mean geometric chord to the aircraft centerline; Q is the dynamic pressure and Sw is the wing area.Then (10) can be integrated into the 6-DOF 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 safety-related parameters. Here, take stall angle of attack (AOA) as an example and the model [20] is shown as (11)ΔCL=CLice-CLαstall=fΔCLHere, αstall 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 high-precision 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 Figure4.Figure 4
Flow chart of flight parameter extremum extraction based on the Monte Carlo simulation.(1) Simulation Status Setting Module. Set the initial flight state of the aircraft, the total number of simulation (n), icing time, and the typical distribution of aircraft icing.(2) Random Variable Extraction Module. The external environment and pilot maneuvering parameters required for the ith 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 Runge-Kutta algorithm. The simulation step length is 0.02s. The flight dynamics model and flight control system are loaded on the real-time 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 ith simulation.(5) Simulation Process Control Module. Let i=i+1, and determine whether i is greater than the setting number n. If not, the simulation returns to the second step. Experiment shows that when n is greater than 2000, the statistical characteristic is stable. So, set n 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 the21st simulation are shown in Figure 5.Figure 5
Flight safety critical parameters variations in the21st simulation.
(a) (b) (c) (d) (e) (f)From Figure5, 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 take-off 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 pilot-in-the-loop 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 Quantiles-Quantiles (QQ) plots, R-square method, correlation coefficient method, and Kolmogorov-Smirnov (K-S) test are performed as Figure 7 and Table 1 show.Table 1
Data goodness of fit for simulation of pilot model and real pilot.
Extreme parameters Correlation coefficient K-S P value R-square Extreme minimum airspeed 0.9607 0.0850 0.4486 0.9230 Extreme roll angle 0.9797 0.0450 0.9855 0.9599 Extreme AOA 0.9667 0.1050 0.2078 0.9346Figure 6
Scatter plot of extremum parameters.Figure 7
QQ diagrams of simulation and test data.The QQ graphical method plots the quantile scatter of the two-sample empirical distribution function on the same graph. The straight lines in Figure7 indicate that the two samples belong to the similar distribution. K-S test is to determine the goodness of fit of the two samples by calculating the empirical probability difference. The P value of the K-S 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 R-square between the samples is greater than 0.9, indicating that the two-sample 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.


## 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 Figure4.Figure 4
Flow chart of flight parameter extremum extraction based on the Monte Carlo simulation.(1) Simulation Status Setting Module. Set the initial flight state of the aircraft, the total number of simulation (n), icing time, and the typical distribution of aircraft icing.(2) Random Variable Extraction Module. The external environment and pilot maneuvering parameters required for the ith 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 Runge-Kutta algorithm. The simulation step length is 0.02s. The flight dynamics model and flight control system are loaded on the real-time 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 ith simulation.(5) Simulation Process Control Module. Let i=i+1, and determine whether i is greater than the setting number n. If not, the simulation returns to the second step. Experiment shows that when n is greater than 2000, the statistical characteristic is stable. So, set n 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 the21st simulation are shown in Figure 5.Figure 5
Flight safety critical parameters variations in the21st simulation.
(a) (b) (c) (d) (e) (f)From Figure5, 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 take-off 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 pilot-in-the-loop 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 Quantiles-Quantiles (QQ) plots, R-square method, correlation coefficient method, and Kolmogorov-Smirnov (K-S) test are performed as Figure 7 and Table 1 show.Table 1
Data goodness of fit for simulation of pilot model and real pilot.
Extreme parameters Correlation coefficient K-S P value R-square Extreme minimum airspeed 0.9607 0.0850 0.4486 0.9230 Extreme roll angle 0.9797 0.0450 0.9855 0.9599 Extreme AOA 0.9667 0.1050 0.2078 0.9346Figure 6
Scatter plot of extremum parameters.Figure 7
QQ diagrams of simulation and test data.The QQ graphical method plots the quantile scatter of the two-sample empirical distribution function on the same graph. The straight lines in Figure7 indicate that the two samples belong to the similar distribution. K-S test is to determine the goodness of fit of the two samples by calculating the empirical probability difference. The P value of the K-S 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 R-square between the samples is greater than 0.9, indicating that the two-sample 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 One-Dimensional 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 Table2, in order to study the distribution characteristics of the flight extreme parameters.Table 2
Statistical characteristics of extreme value samples.
Extreme parameters Maximum Minimum Average Median Variance Kurtosis Skewness Extreme minimum airspeed 115.0038 60.6187 98.9172 102.7124 151.5753 3.2448 -0.9658 Extreme roll angle 103.0343 18.6904 47.4460 44.9050 252.5812 4.0731 0.9875 Extreme AOA 20.9960 4.4056 8.0537 7.1050 8.8444 6.1733 1.7072In Table2, 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 heavy-tailed distribution is common in low-frequency and high risk events. For this type of heavy-tailed 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 heavy-tailed distribution. Lognormal (Logn) distribution, Weibull distribution, and exponential (Exp) distribution can also describe the heavy-tailed characteristics. Here, five typical distribution families are selected to identify one-dimensional extreme value parameters and the normal distribution is also listed to compare, as shown in the following.(12)EV:Fx;μ,σ=-exp⁡-expx-μσ(13)GEV:Fx;ξ,μ,σ=exp⁡-1+ξ·x-μσ-1/ξ.(14)Lognormal:Fx;μ,σ=1σ2π∫0xexp⁡-ln⁡t-μ2/2σ2tdt(15)Weibull:Fx;a,b=∫0xba-btb-1exp⁡-tabdt(16)Exponential:Fx;μ=1-exp⁡-xμ(17)Normal:Fx;μ,σ=1σ2π∫-∞xexp⁡-t-μ22σ2dtThe identified weight parameters for the above equations are listed in Table 3.Table 3
Identified weight parameters.
Model Extreme Vmin Extreme roll angle Extreme AOA EV μ = 104.4189 ,   σ = 8.9329 μ = 56.0386 ,   σ = 19.1285 μ = 9.7361 ,   σ = 4.0631 GEV ξ = - 0.7757 ,   μ = 97.6423, ξ = 0.0012 ,   μ = 40.2577, ξ = 0.3546 ,   μ = 6.5188, σ = 13.5678 σ = 12.3873 σ = 1.4989 Logn μ = 4.5857 ,   σ = 0.1350 μ = 3.8070 ,   σ = 0.3241 μ = 2.0315 ,   σ = 0.3163 Weibull a = 103.9547 ,   b = 10.9225 a = 52.9810 ,   b = 3.0961 a = 9.0375 ,   b = 2.7443 Exp μ = 98.9172 μ = 47.4460 μ = 8.0537 Norm μ = 98.9172 ,   σ = 12.3116 μ = 47.4460 ,   σ = 15.8928 μ = 8.0537 ,   σ = 2.9740In order to judge the accuracy of the above models, the goodness of fit test is applied by the K-S test, Anderson-Darling (A-D) test, andχ2 test. The results are shown in Tables 4–6.Table 4
Fitting test of extreme Vmin.
Model K-S P ( K - S ) A-D P ( A - D ) χ2 EV 0.0831 0.1193 1.8823 0.1068 25.1681 GEV 0.1033 0.0259 2.8973 0.0309 1.8030 Long 0.1465 0.0003 7.5795 0.0002 3.8987 Weibull 0.0426 0.8455 0.3988 0.8498 0.2305 Exp 0.4823 <0.0001 71.1810 <0.0001 25.9954 Norm 0.0831 0.1193 1.8823 0.1068 25.1681Table 5
Fitting test of extreme roll angle.
Model K-S P K - S A-D P ( A - D ) χ2 EV 0.1590 0.0001 10.3930 0 15.1128 GEV 0.0295 0.9931 0.1976 0.9911 0.0921 Long 0.0345 0.9650 0.2201 0.9839 0.0876 Weibull 0.0876 0.0876 3.0570 0.0256 1.4002 Exp 0.3942 <0.0001 43.5374 <0.0001 17.6790 Norm 0.0788 0.1583 2.8350 0.0333 1.2308Table 6
Fitting test of extreme AOA.
Model K-S P K - S A-D P ( A - D ) χ2 EV 0.2606 <0.0001 18.9411 <0.0001 15.5295 GEV 0.0516 0.6413 0.5574 0.6889 0.2089 Long 0.1427 0.0005 4.5021 0.0050 2.4380 Weibull 0.1711 <0.0001 9.8211 <0.0001 5.2385 Exp 0.4548 <0.0001 44.8035 <0.0001 18.1063 Norm 0.1941 <0.0001 10.3963 <0.0001 5.6621In Table4, K-S 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, A-D test is more sensitive to heavy-tailed 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 One-Dimensional 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αmax/αcδf,Ma,ice. 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 (18)P=1,ifVminVcδf,ice<1P=1,ifϕmax85>1P=1,ifαmaxαcδf,Ma,ice>1(19)Pv=F1V¯min<1(20)Pϕ=1-F2ϕ¯max>1(21)Pα=1-F3α¯max>1Based 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 MIL-STD-882D, 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 three-dimensional flight extremum parameters with different mathematical distributions is built.


## 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 Table2, in order to study the distribution characteristics of the flight extreme parameters.Table 2
Statistical characteristics of extreme value samples.
Extreme parameters Maximum Minimum Average Median Variance Kurtosis Skewness Extreme minimum airspeed 115.0038 60.6187 98.9172 102.7124 151.5753 3.2448 -0.9658 Extreme roll angle 103.0343 18.6904 47.4460 44.9050 252.5812 4.0731 0.9875 Extreme AOA 20.9960 4.4056 8.0537 7.1050 8.8444 6.1733 1.7072In Table2, 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 heavy-tailed distribution is common in low-frequency and high risk events. For this type of heavy-tailed 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 heavy-tailed distribution. Lognormal (Logn) distribution, Weibull distribution, and exponential (Exp) distribution can also describe the heavy-tailed characteristics. Here, five typical distribution families are selected to identify one-dimensional extreme value parameters and the normal distribution is also listed to compare, as shown in the following.(12)EV:Fx;μ,σ=-exp⁡-expx-μσ(13)GEV:Fx;ξ,μ,σ=exp⁡-1+ξ·x-μσ-1/ξ.(14)Lognormal:Fx;μ,σ=1σ2π∫0xexp⁡-ln⁡t-μ2/2σ2tdt(15)Weibull:Fx;a,b=∫0xba-btb-1exp⁡-tabdt(16)Exponential:Fx;μ=1-exp⁡-xμ(17)Normal:Fx;μ,σ=1σ2π∫-∞xexp⁡-t-μ22σ2dtThe identified weight parameters for the above equations are listed in Table 3.Table 3
Identified weight parameters.
Model Extreme Vmin Extreme roll angle Extreme AOA EV μ = 104.4189 ,   σ = 8.9329 μ = 56.0386 ,   σ = 19.1285 μ = 9.7361 ,   σ = 4.0631 GEV ξ = - 0.7757 ,   μ = 97.6423, ξ = 0.0012 ,   μ = 40.2577, ξ = 0.3546 ,   μ = 6.5188, σ = 13.5678 σ = 12.3873 σ = 1.4989 Logn μ = 4.5857 ,   σ = 0.1350 μ = 3.8070 ,   σ = 0.3241 μ = 2.0315 ,   σ = 0.3163 Weibull a = 103.9547 ,   b = 10.9225 a = 52.9810 ,   b = 3.0961 a = 9.0375 ,   b = 2.7443 Exp μ = 98.9172 μ = 47.4460 μ = 8.0537 Norm μ = 98.9172 ,   σ = 12.3116 μ = 47.4460 ,   σ = 15.8928 μ = 8.0537 ,   σ = 2.9740In order to judge the accuracy of the above models, the goodness of fit test is applied by the K-S test, Anderson-Darling (A-D) test, andχ2 test. The results are shown in Tables 4–6.Table 4
Fitting test of extreme Vmin.
Model K-S P ( K - S ) A-D P ( A - D ) χ2 EV 0.0831 0.1193 1.8823 0.1068 25.1681 GEV 0.1033 0.0259 2.8973 0.0309 1.8030 Long 0.1465 0.0003 7.5795 0.0002 3.8987 Weibull 0.0426 0.8455 0.3988 0.8498 0.2305 Exp 0.4823 <0.0001 71.1810 <0.0001 25.9954 Norm 0.0831 0.1193 1.8823 0.1068 25.1681Table 5
Fitting test of extreme roll angle.
Model K-S P K - S A-D P ( A - D ) χ2 EV 0.1590 0.0001 10.3930 0 15.1128 GEV 0.0295 0.9931 0.1976 0.9911 0.0921 Long 0.0345 0.9650 0.2201 0.9839 0.0876 Weibull 0.0876 0.0876 3.0570 0.0256 1.4002 Exp 0.3942 <0.0001 43.5374 <0.0001 17.6790 Norm 0.0788 0.1583 2.8350 0.0333 1.2308Table 6
Fitting test of extreme AOA.
Model K-S P K - S A-D P ( A - D ) χ2 EV 0.2606 <0.0001 18.9411 <0.0001 15.5295 GEV 0.0516 0.6413 0.5574 0.6889 0.2089 Long 0.1427 0.0005 4.5021 0.0050 2.4380 Weibull 0.1711 <0.0001 9.8211 <0.0001 5.2385 Exp 0.4548 <0.0001 44.8035 <0.0001 18.1063 Norm 0.1941 <0.0001 10.3963 <0.0001 5.6621In Table4, K-S 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, A-D test is more sensitive to heavy-tailed 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 One-Dimensional 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αmax/αcδf,Ma,ice. 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 (18)P=1,ifVminVcδf,ice<1P=1,ifϕmax85>1P=1,ifαmaxαcδf,Ma,ice>1(19)Pv=F1V¯min<1(20)Pϕ=1-F2ϕ¯max>1(21)Pα=1-F3α¯max>1Based 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 MIL-STD-882D, 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 three-dimensional flight extremum parameters with different mathematical distributions is built.


## 5. Three-Dimensional 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 vectorVmin,ϕmax,αmax is FVmin,ϕmax,αmax and the marginal distribution functions are F1Vmin, F2ϕmax, F3αmax. We have concluded that F1Vmin obeys Weibull distribution and F2ϕmax and F3αmax obey GEV distribution at 95% confidence level in Section 4. According to Sklar theorem [8], for any Vmin,ϕmax,αmax∈Rd, there must be a Copula function C, which satisfies (22)CF1Vmin,F2ϕmax,F3αmax=FVmin,ϕmax,αmaxA three-dimensional Copula model is constructed based on asymmetric Archimedes Copula, as shown in(23)Cu,v,w=C1w,C2u,v=φ1-1φ1w+φ1φ2-1φ2u+φ2vwhere φ(·) is Archimedes Copula generating function; u,v,w is marginal distribution of the sample Vmin,ϕmax,αmax. Firstly, Copula C2 is generated by u and v; secondly, Copula C1 is generated by combining w and C2. Thirdly, assuming θ1 is the parameter of Copula C2 and θ2 is the parameter of Copula C1, expressions for C2(u,v;θ1) and C1(w,C2;θ2) are given; finally, according to (23), C2(u,v;θ1) is brought into C1(w,C2;θ2). Then, the three-dimensional 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 safetyHere,F1Vmin, F2ϕmax, F3αmax are all continuous distribution functions, so Copula C is the only one. The commonly used multidimensional Copula models and corresponding generators are shown as (24)-(28).Gumbel Copula:(24)Cu,v,w=exp⁡--ln⁡wθ2+-ln⁡uθ1+-ln⁡vθ1θ2/θ11/θ2φt=-ln⁡tθFrank Copula:(25)Cu,v,w=-1θ2log⁡1-1-e-θ2-11-e-θ2w1-1-1-e-θ1-11-e-θ1u1-e-θ1vθ2/θ1φt=ln1-e-θt1-e-θClayton Copula:(26)Cu,v,w=w-θ2+u-θ1+v-θ1-1θ2/θ1-1-1/θ2φt=1θt-θ-1GS Copula:(27)Cu,v,w=1+1w-1θ2+1u-1θ1+1v-1θ1θ2/θ11/θ2-1φt=1t-1θJoe Copula:(28)Cu,v,w=1-1-uθ11-1-vθ1+1-vθ1θ2/θ11-1-wθ2+1-wθ21/θ2φt=-ln⁡1-1-tθ


### 5.2. Identification of the Parameters in Copula Functions

Using nonlinear identification algorithms, the identified weight parametersθ1 and θ2 for the above models are calculated and listed in Table 7.Table 7
Identified weight parameters.
Copula model Weight parameters Clayton Copula θ 1 = 9.2620 θ 2 = 9.0787 Gumbel Copula θ 1 = 1.7251 θ 2 = 2.4825 Frank Copula θ 1 = 40.0621 θ 2 = 35.3235 Gs Copula θ 1 = 1.5518 θ 2 = 1.9015 Joe Copula θ 1 = 12.6297 θ 2 = 10.6006Goodness of fit test is applied by the Akaike Information Criteria (AIC), Bayesian Information Criterion (BIC),χ2 test, and K-S test, in order to judge the Copula model accuracy. The results are shown in Table 8. We can analyze that the PK-S 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 chi-square statistics, Joe model can more accurately describe the three-dimensional extreme value distribution model under the condition than Clayton model. What is more, Joe model is sensitive to heavy-tail characteristic. So we choose Joe model to describe the joint probability distribution of minimum airspeed extreme, maximum roll angle extreme, and maximum AOA extreme.Table 8
Fitting test for different Copula models.
Copula model AIC BIC χ2 K-S P K - S Clayton Copula -1124.1 -1117.5 8.2103 0.1100 0.1668 Gumbel Copula -771.97 -765.37 45035 0.6500 4.3188E-38 Frank Copula -9.6640 -6.6534 0.4119 0.2986 0.0204 Gs Copula -631.4853 -624.887 11.0808 0.2000 5.4431E-4 Joe Copula -1416.8 -1410.2 1.9402 0.1050 0.2078


### 5.3. Three-Dimensional Extreme Parameter Flight Risk Probability and Analysis

According to the risk discriminant (18) and Joe Copula model, the flight risk probability of three-dimensional extreme parameters is shown in (29)Pr=1-CJoe(F1VminVcδf,ice<1,F2ϕmax85>1,F3αmaxαcδf,Ma,ice>1.C J o e 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 one-dimensional 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 one-dimensional 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, Figure8 shows the joint probability density distribution of Joe Copula model and Clayton Copula model when w=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 cliff-breaking 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.Figure 8
Probability density maps of Joe and Clayton copula models (w=0.8).


## 5.1. Mathematical Distribution Type Hypothesis

Assume the distribution function of the extremum vectorVmin,ϕmax,αmax is FVmin,ϕmax,αmax and the marginal distribution functions are F1Vmin, F2ϕmax, F3αmax. We have concluded that F1Vmin obeys Weibull distribution and F2ϕmax and F3αmax obey GEV distribution at 95% confidence level in Section 4. According to Sklar theorem [8], for any Vmin,ϕmax,αmax∈Rd, there must be a Copula function C, which satisfies (22)CF1Vmin,F2ϕmax,F3αmax=FVmin,ϕmax,αmaxA three-dimensional Copula model is constructed based on asymmetric Archimedes Copula, as shown in(23)Cu,v,w=C1w,C2u,v=φ1-1φ1w+φ1φ2-1φ2u+φ2vwhere φ(·) is Archimedes Copula generating function; u,v,w is marginal distribution of the sample Vmin,ϕmax,αmax. Firstly, Copula C2 is generated by u and v; secondly, Copula C1 is generated by combining w and C2. Thirdly, assuming θ1 is the parameter of Copula C2 and θ2 is the parameter of Copula C1, expressions for C2(u,v;θ1) and C1(w,C2;θ2) are given; finally, according to (23), C2(u,v;θ1) is brought into C1(w,C2;θ2). Then, the three-dimensional 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 safetyHere,F1Vmin, F2ϕmax, F3αmax are all continuous distribution functions, so Copula C is the only one. The commonly used multidimensional Copula models and corresponding generators are shown as (24)-(28).Gumbel Copula:(24)Cu,v,w=exp⁡--ln⁡wθ2+-ln⁡uθ1+-ln⁡vθ1θ2/θ11/θ2φt=-ln⁡tθFrank Copula:(25)Cu,v,w=-1θ2log⁡1-1-e-θ2-11-e-θ2w1-1-1-e-θ1-11-e-θ1u1-e-θ1vθ2/θ1φt=ln1-e-θt1-e-θClayton Copula:(26)Cu,v,w=w-θ2+u-θ1+v-θ1-1θ2/θ1-1-1/θ2φt=1θt-θ-1GS Copula:(27)Cu,v,w=1+1w-1θ2+1u-1θ1+1v-1θ1θ2/θ11/θ2-1φt=1t-1θJoe Copula:(28)Cu,v,w=1-1-uθ11-1-vθ1+1-vθ1θ2/θ11-1-wθ2+1-wθ21/θ2φt=-ln⁡1-1-tθ


## 5.2. Identification of the Parameters in Copula Functions

Using nonlinear identification algorithms, the identified weight parametersθ1 and θ2 for the above models are calculated and listed in Table 7.Table 7
Identified weight parameters.
Copula model Weight parameters Clayton Copula θ 1 = 9.2620 θ 2 = 9.0787 Gumbel Copula θ 1 = 1.7251 θ 2 = 2.4825 Frank Copula θ 1 = 40.0621 θ 2 = 35.3235 Gs Copula θ 1 = 1.5518 θ 2 = 1.9015 Joe Copula θ 1 = 12.6297 θ 2 = 10.6006Goodness of fit test is applied by the Akaike Information Criteria (AIC), Bayesian Information Criterion (BIC),χ2 test, and K-S test, in order to judge the Copula model accuracy. The results are shown in Table 8. We can analyze that the PK-S 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 chi-square statistics, Joe model can more accurately describe the three-dimensional extreme value distribution model under the condition than Clayton model. What is more, Joe model is sensitive to heavy-tail characteristic. So we choose Joe model to describe the joint probability distribution of minimum airspeed extreme, maximum roll angle extreme, and maximum AOA extreme.Table 8
Fitting test for different Copula models.
Copula model AIC BIC χ2 K-S P K - S Clayton Copula -1124.1 -1117.5 8.2103 0.1100 0.1668 Gumbel Copula -771.97 -765.37 45035 0.6500 4.3188E-38 Frank Copula -9.6640 -6.6534 0.4119 0.2986 0.0204 Gs Copula -631.4853 -624.887 11.0808 0.2000 5.4431E-4 Joe Copula -1416.8 -1410.2 1.9402 0.1050 0.2078


## 5.3. Three-Dimensional Extreme Parameter Flight Risk Probability and Analysis

According to the risk discriminant (18) and Joe Copula model, the flight risk probability of three-dimensional extreme parameters is shown in (29)Pr=1-CJoe(F1VminVcδf,ice<1,F2ϕmax85>1,F3αmaxαcδf,Ma,ice>1.C J o e 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 one-dimensional 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 one-dimensional 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, Figure8 shows the joint probability density distribution of Joe Copula model and Clayton Copula model when w=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 cliff-breaking 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.Figure 8
Probability density maps of Joe and Clayton copula models (w=0.8).


## 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 real-time 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.



---

*Source: 2546395-2019-04-09.xml*