A civil aircraft operation risk assessment method for multiple failure modes

Abstract

The present invention discloses a method for assessing the operation risk of civil aircraft with multiple failure modes in the field of aircraft safety technology, comprising: aiming at the characteristics of multiple wear and tear failures of civil aircraft structural parts, establishing a failure model based on mixed Weibull distribution, and estimating model parameters based on the EM algorithm optimized by PSO; executing a Monte Carlo simulation algorithm, simulating the number DA of defective aircraft with multiple wear and tear failures in a fleet according to the failure model; constructing a Bayesian network to perform failure propagation analysis on the initial cause events leading to unsafe results, and calculating the probability CP of unsafe results under the condition of the initial cause events. The present invention carries out operation risk assessment for multiple failure modes, overcoming the defect that the existing risk assessment method can only be applied to a single failure mode.

Description

Technical Field

[0001] The present invention relates to the field of aircraft safety technology, and in particular to a civil aircraft operation risk assessment method oriented to multiple failure modes. Background Art

[0002] Civilian aircraft have become the preferred means of transportation for more and more people due to their high efficiency. For the civil aviation transportation industry, safety is the top priority. It is very important for all civil aircraft in service to ensure that aircraft can continue to fly safely under the expected operating environment and usage restrictions.

[0003] The existing aviation safety analysis methods are mainly aimed at individual aircraft, and a relatively complete safety management system including quantitative and qualitative analysis has been formed. The Transport Airplane Risk Assessment Methodology (TARAM) is a quantitative risk analysis method in the Federal Aviation Administration (FAA) Directive 8110.107A Monitor Safety-Analyze Data (MASD). This method provides a quantitative method for different operating periods and different risk values ​​for fleets and single aircraft. It is used to implement the quantitative risk analysis of the continuous operational safety of the transport fleet required by the MSAD directive. For different periods of fleet operation, the method provides three types of wear and tear failure risk values: overall uncorrected fleet risk, 90-day fleet risk, and control procedure fleet risk. Among them, the overall uncorrected fleet risk (R T ) is used to make long-term predictions on the risk of unsafe conditions in the fleet and guide safety decisions. R can be calculated according to the formula T =DA×ND×CP×IR. IR and ND are usually given by fleet operation experience. For CP, it is necessary to establish a causal chain from the initial event to the unsafe result for analysis, and combine historical data to obtain the probability of each failure event. The final probability value is obtained by accumulating the conditional probabilities of all conditions of the unsafe result. For the case where the failure mode of the initial event is unknown, the Weibull analysis tool or reliability analysis software can be used to obtain the failure distribution parameters. According to the formula To calculate the DA value, the existing civil aircraft operation risk assessment method is to study the failure rate or single wear and tear failure by using the TARAM method.

[0004] However, the inventors of this application found that the above technology has at least the following technical problems in the process of implementing the technical solution of the invention in the embodiment of this application:

[0005] 1. The civil aircraft risk analysis method only considers the impact of a single failure on aircraft operation safety. However, as the structure of civil aircraft becomes more complex and the service time increases, the wear and tear failures usually present the characteristics of multiple failure modes. This method cannot accurately assess the fleet operation risk.

[0006] 2. When multiple failure modes occur in the fleet, the corresponding failure distribution function is relatively complex. The analytical method needs to calculate the failure cumulative distribution function value of each aircraft's current time and retirement time multiple times. The calculation process is complicated and inefficient.

[0007] 3. When studying the failure propagation relationship from the initial cause event to the unsafe result, the risk value of each branch is calculated one by one through causal chain analysis, and the final CP value is accumulated. When there are too many branches in the causal chain, the calculation process is extremely cumbersome.

[0008] Based on this, the present invention designs a civil aircraft operation risk assessment method for multiple failure modes to solve the above problems. Summary of the invention

[0009] The purpose of the present invention is to provide a civil aircraft operation risk assessment method for multiple failure modes to solve the technical problems mentioned in the existing background technology proposed in the above background technology.

[0010] To achieve the above object, the present invention provides the following technical solution: a method for assessing the operation risk of a civil aircraft for multiple failure modes, comprising:

[0011] Considering the characteristics of multiple wear and failure of civil aircraft structural parts, a failure model of aircraft structural parts based on mixed Weibull distribution is established by using the failure data of fleet operation, and the model parameters are estimated based on the EM algorithm optimized by PSO.

[0012] Execute Monte Carlo simulation algorithm to simulate the number of defective aircraft with multiple wear and tear failures in the fleet according to the failure model;

[0013] Construct a Bayesian network to perform failure propagation analysis on the initial cause events that lead to unsafe results, and calculate the probability CP of unsafe results under the condition of the initial cause events;

[0014] According to the predetermined damage rate IR and non-detection rate ND, the fleet risk level R is calculated. T =DA×ND×CP×IR.

[0015] Furthermore, the construction of the aircraft structural component failure model includes:

[0016] Using failure data T = [t i ], construct the likelihood function of complete data:

[0017]

[0018] Where N is the total amount of sample data, K is the number of Weibull sub-distributions, and Z = [Z ik ],Z ik is a hidden variable, and when the i-th observation comes from the k-th Weibull subdistribution, then Z ik =1, otherwise Z ik =0; θ is the parameter set of the distribution, expressed as θ=(π1, ..., π K ; θ1, ..., θ K ), π k is the kth sub-distribution weighting factor, θ k =(α k , β k ) represents the scale parameter α of the kth sub-distribution k and shape parameter β k ;f(·|θ k ) is the probability density function of the kth Weibull subdistribution;

[0019] Execute the EM algorithm by giving the failure data T and the current parameter estimate θ (j) Find the conditional probability P(Z|T,θ) of Z (j) ), and obtain the expected value Q of the likelihood function of the conditional probability:

[0020]

[0021] Taking the Q function as the objective function, the optimization model for mixed Weibull distribution parameter estimation is established as follows:

[0022] maxQ(θ,θ (j) )

[0023]

[0024] Furthermore, the process of executing the EM algorithm also includes:

[0025] Determine the parameter θ (j) The Q function corresponding to the jth iteration value and Q(θ,θ (j) ) is maximized to obtain the j+1th iteration result θ (j+1) , repeat the above steps until convergence.

[0026] Furthermore, the Monte Carlo simulation algorithm is executed to simulate the number of defective aircraft with multiple wear and tear failures in the fleet according to the failure model. The number DA includes:

[0027] Calculate the cumulative probability of the weighting factor Where q0 = 0, k is the number of Weibull subdistributions;

[0028] Generate a random number r from a uniform distribution between [0, 1] for each aircraft and determine the probability interval [q k-1 ,q k ], and the aircraft is found to obey the kth Weibull sub-distribution failure mode;

[0029] Perform inverse transformation simulation to generate aircraft failure time Among them, α k , β k are the scale parameter of the kth sub-distribution and the shape parameter of the kth sub-distribution respectively;

[0030] Statistics of all aircraft in the fleet during the failure prediction period T m quantity, and T m The average value after multiple simulations is taken as the DA value of the number of defective aircraft.

[0031] Furthermore, the Bayesian network is constructed to perform failure propagation analysis on the initial events that lead to unsafe results, including:

[0032] Analyze the failure propagation chain of aircraft structural parts from initial cause events to unsafe consequences;

[0033] Each level of intermediate failure corresponds to each node event in the BN;

[0034] Specific symbols are selected to describe the causal relationship between nodes, and a directed acyclic graph is constructed to intuitively describe the development of unsafe results.

[0035] Furthermore, the calculation of the probability of occurrence of unsafe results CP under the condition of the occurrence of the initial cause event includes:

[0036] The BN joint distribution is decomposed into the product form of factor functions through variable elimination method;

[0037] Obtain a factor product formula containing only the query variable by eliminating the corresponding variables in each factor function in an elimination order;

[0038] The initial event is set as the evidence variable, and the unsafe result is set as the query variable. The distribution of each node is updated through variable elimination and the CP value of the unsafe result under the condition of the initial event occurring is calculated.

[0039] Furthermore, the optimization of the EM algorithm by PSO includes:

[0040] Correspond the D-dimensional position information of each particle in PSO to the parameter θ;

[0041] When the algorithm is initialized, the update limits of particle position and velocity are determined according to the initial values ​​of each parameter, and a particle population of popsize is randomly generated;

[0042] The Q function is used as the fitness function of PSO. The particles are updated by comparing the fitness values ​​of each generation of particles, and the particles with the best fitness values ​​are retained.

[0043] Update the particles to the maximum number of optimization times and output the global optimal particle position as the maximization result of this round of EM iteration.

[0044] One or more technical solutions provided in the embodiments of the present invention have at least the following technical effects or advantages:

[0045] 1. Carry out operational risk assessment for multiple failure modes, overcoming the defect that existing risk assessment methods can only be applied to a single failure mode;

[0046] 2. By using BN to analyze the failure propagation relationship of the initial cause event, the tedious process of using causal chain analysis when multiple intermediate failures occur simultaneously is avoided;

[0047] 3. The EM algorithm based on PSO optimization overcomes the shortcomings of commonly used search and approximation algorithms that are difficult to solve the problem of mixed distribution parameter estimation. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings required for describing the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other accompanying drawings can be obtained based on these accompanying drawings without paying creative work.

[0049] Figure 1 is a flowchart of the evaluation method of the present invention;

[0050] Figure 2 It is a flowchart of the EM algorithm based on PSO optimization of the present invention;

[0051] Figure 3 A probability density diagram of failure of first-cause events of a fleet implementing the present invention;

[0052] Figure 4 A current operating status diagram of the fleet for implementing the present invention;

[0053] Figure 5 This is a flowchart of the DA process based on Monte Carlo simulation of the fleet in the present invention;

[0054] Figure 6 This is a diagram of the DA simulation results of the fleet implemented in the present invention;

[0055] Figure 7 BN diagrams before and after the simplification of the wing cracks of the fleet according to the present invention;

[0056] Figure 8A data table showing the life of wing ribs of a fleet implementing the present invention;

[0057] Fig. 9 This is a diagram of the fleet BN node event information table implemented in the present invention;

[0058] Fig.10 This is a table showing damage rates of different results for a fleet implementing the present invention. DETAILED DESCRIPTION

[0059] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0060] See also Figure 1-10 The embodiment of the present invention provides a technical solution: a method for assessing the operation risk of a civil aircraft for multiple failure modes, comprising:

[0061] Considering the characteristics of multiple wear and failure of civil aircraft structural parts, a failure model of aircraft structural parts based on mixed Weibull distribution is established by using the failure data of fleet operation, and the model parameters are estimated based on the EM algorithm optimized by PSO.

[0062] Execute Monte Carlo simulation algorithm to simulate the number of defective aircraft with multiple wear and tear failures in the fleet according to the failure model;

[0063] Construct a Bayesian network to perform failure propagation analysis on the initial cause events that lead to unsafe results, and calculate the probability CP of unsafe results under the condition of the initial cause events;

[0064] According to the predetermined damage rate IR and non-detection rate ND, the fleet risk level R is calculated. T =DA×ND×CP×IR.

[0065] Through the above steps, it is not difficult to see that in the civil aircraft operation risk assessment method for multiple failure modes of the present application, the failure data in the fleet operation is used to establish a multiple failure model of aircraft structural parts based on the mixed Weibull distribution, and based on the PSO optimization EM algorithm, parameter estimation of the multiple failure model can be achieved. Moreover, based on the PSO optimization EM algorithm, it can effectively overcome the shortcomings of commonly used search and approximation algorithms that are difficult to solve the mixed distribution parameter estimation problem.

[0066] Weibull distribution is the most widely used life distribution in reliability models. Its main advantages are good data fitting, suitability for small sample sampling, and strong adaptability to various types of data. The mixed Weibull distribution is used to fit the life distribution of failure data, which can realize the construction of the likelihood function of complete data. The construction of the aircraft structural failure model includes:

[0067] Using failure data T = [t i ], construct the likelihood function of complete data:

[0068]

[0069] Where N is the total amount of sample data, K is the number of Weibull sub-distributions, and Z = [Z ik ],Z ik is a hidden variable, and when the i-th observation comes from the k-th Weibull subdistribution, then Z ik =1, otherwise Z ik =0; θ is the parameter set of the distribution, expressed as θ=(π1, ..., π K ; θ1, ..., θ K ), π k is the kth sub-distribution weighting factor, θ k =(α k , β k ) represents the scale parameter α of the kth sub-distribution k and shape parameter β k ;f(·|θ k ) is the probability density function of the kth Weibull subdistribution;

[0070] Execute the EM algorithm by giving the failure data T and the current parameter estimate θ (j) Find the conditional probability P(Z|T,θ) of Z (j) ), and obtain the expected value Q of the likelihood function of the conditional probability:

[0071]

[0072] Taking the Q function as the objective function, the optimization model for mixed Weibull distribution parameter estimation is established as follows:

[0073] maxQ(θ,θ (j) )

[0074]

[0075] Since the Q function of the mixed Weibull distribution parameters is usually a complex nonlinear equation, its analytical value cannot be directly obtained. When the EM algorithm based on PSO optimization is used to implement the process of mixed Weibull distribution parameter estimation, the execution of the EM algorithm also includes:

[0076] Determine the parameter θ (j) The Q function corresponding to the jth iteration value and Q(θ,θ (j) ) is maximized to obtain the j+1th iteration result θ (j+1) , repeat the above operation until convergence;

[0077] It should be noted that EM is a maximum likelihood estimation method for probability distribution parameters with hidden variables, which mainly includes: expectation step (E-Step) and maximization step (M-Step). In the E-step, the j-th iteration parameter result θ is determined (j) The corresponding Q function value, and in the M step Q(θ, θ (j) ) is maximized to obtain the j+1th iteration result θ (j+1) , repeat the above two steps until the convergence condition is reached.

[0078] Furthermore, the Monte Carlo simulation algorithm is executed to simulate the number of defective aircraft with multiple wear and tear failures in the fleet through the failure model. The number DA includes:

[0079] Calculate the cumulative probability of the weighting factor Where q0 = 0, k is the number of Weibull subdistributions;

[0080] Generate a random number r from a uniform distribution between [0, 1] for each aircraft and determine the probability interval [q k-1 ,q k ], and the aircraft is found to obey the kth Weibull sub-distribution failure mode;

[0081] Perform inverse transformation simulation to generate aircraft failure time where α k , β k are the scale parameter of the kth sub-distribution and the shape parameter of the kth sub-distribution respectively;

[0082] Statistics of all aircraft in the fleet during the failure prediction period T m quantity, and T m The average value after multiple simulations is taken as the DA value of defective aircraft quantity;

[0083] When executing the Monte Carlo simulation algorithm to simulate the number of defective aircraft DA with multiple wear and tear failures in the fleet through the failure model, 1. Input the total number of aircraft in the fleet M and the current operating time T of each aircraft. 0(mm) and retirement time T R , determine the maximum number of loops iter max , initialize DA and start the loop; 2. Input the parameter estimation results of the multiple failure distribution model, and calculate the cumulative probability of each sub-distribution weighting factor for the Weibull sub-distribution with different weights Where q0 = 0, randomly generate a random number r that follows a uniform distribution between [0, 1], and calculate the probability interval [q k-1 ,q k ], determine that the mth aircraft obeys the kth Weibull sub-distribution failure mode, and simulate the failure time of the aircraft to be 3. Determine the aircraft T m Is it during the failure prediction period [T 0(m) , T R ], if yes, then the DA number is increased by 1, otherwise return to step 2 to generate the next aircraft failure time; 4. Traverse each aircraft in the fleet to complete a simulation and record the DA value obtained in the simulation; repeat 2 to 4 until the maximum number of loops iter max , output the average value of all DA.

[0084] The BN implementation method for conducting failure propagation analysis on the initial cause event leading to the unsafe result is to construct a Bayesian network to conduct failure propagation analysis on the initial cause event leading to the unsafe result, including:

[0085] Analyze aircraft structural components and failure propagation chains leading to unsafe outcomes;

[0086] Each level of intermediate failure corresponds to each node event in the BN;

[0087] Select specific symbols to describe the causal relationship between nodes and construct a directed acyclic graph that intuitively describes the development of unsafe results;

[0088] It should be noted that, through the above content, the failure propagation relationship of aircraft structural parts from the initial cause event to the final unsafe result is first analyzed. Each level of intermediate failure corresponds to the node event in the BN, and the causal relationship between the nodes is represented by "edges", so as to construct a directed acyclic graph that intuitively describes the development of unsafe results.

[0089] When using the variable elimination method to reason about BN, the probability CP of unsafe results under the condition of the occurrence of the initial cause event is calculated, including:

[0090] The BN joint distribution is decomposed into the product form of factor functions through variable elimination method;

[0091] Obtain a factor product formula containing only the query variable by eliminating the corresponding variables in each factor function in an elimination order;

[0092] The initial event is set as the evidence variable, and the unsafe result is set as the query variable. The distribution of each node is updated through variable elimination, and the CP value of the unsafe result under the condition of the initial event is calculated.

[0093] Specifically, the BN joint distribution is decomposed into the product form of factor functions to obtain a probability distribution set ζ composed of factor functions; the initial event is set as the evidence variable and the unsafe result is set as the query variable, the observation value is set for the evidence variable and the factor formula containing the evidence is updated; the remaining variables except the evidence and query variables are sorted, and the contained variables are marginalized and eliminated in order in the updated ζ; the factor formulas in ζ after elimination are multiplied to obtain the probability distribution containing only the query variable, and then normalized to obtain the posterior probability of the required variable; in risk assessment, the initial event is set as the evidence variable and the unsafe result is set as the query variable, and the variable elimination method is used to update the distribution of each node to obtain the CP value of the final unsafe result.

[0094] Since the Q function of the mixed Weibull distribution parameters is usually a complex nonlinear equation, its analytical value cannot be directly obtained. Therefore, the EM algorithm optimized by PSO includes:

[0095] Correspond the D-dimensional position information of each particle in PSO to the parameter θ;

[0096] When the algorithm is initialized, the update limits of particle position and velocity are determined according to the initial values ​​of each parameter, and a particle population of popsize is randomly generated;

[0097] The Q function is used as the fitness function of PSO, and the particles are updated by comparing the fitness values ​​of each generation of particles;

[0098] Determine whether the updated particle has reached the maximum number of optimization times, and when the maximum number of optimization times is reached, output the global optimal particle position as the maximization result.

[0099] The specific implementation process of this application is:

[0100] like Figure 1 As shown, S1. Combined with the aircraft component failure data obtained during fleet operation, a mixed Weibull distribution parameter estimation model based on the optimized EM algorithm is constructed. Taking a passenger fleet with a total of 500 aircraft as an example, an MSAD report pointed out that the fleet found cracks in the wing structures of many aircraft during regular maintenance. After preliminary analysis, it is highly likely that the cracks in the fleet are caused by fatigue damage due to two reasons. As the service time on each aircraft accumulates, the growth of wing rib cracks to fracture will cause excess loads to be transferred to the surrounding structural components, which will eventually lead to serious consequences such as aircraft loss of control. Therefore, it is necessary to conduct a risk assessment of the affected fleet. From historical experience, it is known that the shape parameter of the fatigue crack failure is β = 2.78, so the constructed double mixed Weibull distribution parameter EM model can be simplified to:

[0101] maxQ(θ,θ (j) )

[0102]

[0103] The life samples of crack growth to failure that are currently available are as follows: Figure 8 As shown in Table 1. The life data is input into the parameter estimation model, using Figure 2 The EM algorithm process based on PSO optimization is given to solve the model. The specific steps are as follows:

[0104] (1) Parameter and particle initialization: The D-dimensional position information of each particle in the PSO corresponds to the required estimated parameter θ. The approximate value of each parameter is determined by the linearization method, and the particle position limit is set according to the initial value of each parameter. and The particle population number popsize and the maximum number of particle optimization cycles maxiter. In this case, the parameter vector θ = (α1, α2, π) is used as the three-dimensional position information of the particle. The particle position and velocity initialization formula of the dth dimension is

[0105] x d =x d +rand·(x d max -x d min )

[0106] v d =0.1·(x d max -x d min );

[0107] The obtained particle position is set as the initial value of the particle individual optimal position pbest and the global optimal position gbest at the same time;

[0108] (2) E-step: In the expected step, the hidden variable Z and Q function need to be determined by the current parameters. First, the position and velocity of each particle are updated. The velocity and position update formula of the sth generation particle is:

[0109]

[0110]

[0111] In the formula, c1 and c2 are adjustment coefficients, whose purpose is to adjust the role of the particle's own experience and global experience in motion. In this case, the acceleration coefficients are both 0.5; r1 and r2 are random numbers between (0, 1); w is the inertia weight, which is commonly linear adjustment weight and random adjustment weight. This model uses random adjustment method.

[0112] w=0.5+rand / 2.0

[0113] Calculate the hidden variable estimate based on the model parameters generated by the particles for:

[0114]

[0115]

[0116] Among them, t i is the input case failure sample, the total number of data is 20. Let θ (j) is the estimated value of the parameter obtained after the jth EM iteration, and is substituted into θ in the next iteration. (j) and Calculate Q(θ,θ (j) )value.

[0117] (3) M-step: The core problem of the maximization step is to solve Q(θ, θ (j) ) to maximize the parameter θ and obtain the estimated value of the model parameter for the next round of iteration. The objective function of PSO optimization is:

[0118]

[0119] By comparing the size of the Q function value of each generation of particles, the individual and global positions are updated. The particles are optimized to the maximum number of cycles. The obtained global optimal position of the particles is the estimated result θ of this round of EM iteration. (j+1) In order to improve the estimation accuracy of the weighting factor, it can be obtained based on the current latent variables:

[0120]

[0121] Compare The size of the Q function value before and after the update is taken as the input weighting factor for the next round of iteration.

[0122] (4) Determine the end condition of the iteration: Repeat steps (2) and (3) until the convergence condition is reached

[0123] ||Q(θ j+1 ,θ (j) )-Q(θ j ,θ (j) )||<ε

[0124] In the formula, ε is a small positive number. After the convergence condition is met, the optimal estimated parameter value of the model is output.

[0125] The estimated values ​​of the distribution scale parameters obtained by the above algorithm are The estimated weighting factor is The probability density function of the crack failure event of the fleet is as follows: Figure 3As shown in the figure, the scattered points represent the probability density values ​​of the true failure sample values. It should be noted that when α is less than the aircraft retirement time, it indicates that the damage occurs in the early stage of life. On the contrary, when α is greater than the retirement time, it indicates that the damage develops slowly.

[0126] S2. Based on the obtained distribution parameters, the Monte Carlo simulation algorithm is used to simulate the number of DAs in the fleet affected by multiple wear and tear failure modes. The specific definition of DA in fleet risk is the weighted events or deaths that are expected to occur within a specified time (usually the time from the current service age of the aircraft to the retirement age) if no measures are taken to correct the identified potential unsafe conditions. In this case, the number of aircraft in the fleet M is 500, and the current service age of each aircraft is as follows: Figure 4 As shown in the figure, the age unit is flight cycle (FC), and the average retirement age of the fleet is T R =50000FC. Input the known parameters into the Monte Carlo simulation algorithm, such as Figure 5 As shown, set iter max The number of cycles is 5000, DA is initialized to 0 and the cycle starts. In this case, the cumulative probabilities of the weighted factors of each sub-distribution are q1=0.26, q1=1, and the Weibull sub-distribution failure mode of the mth aircraft is determined according to the probability interval [q0, q1] or [q1, q2] to which r belongs. The inverse transformation method is used to generate the failure time T of the aircraft structure. m , to determine whether it is within the failure prediction period [T 0(m) , T R ], the final simulation result is shown in the figure Figure 6 As shown, the simulated fleet DA value is approximately 153 aircraft.

[0127] S3. Comprehensively analyze the conditional probability of unsafe results in the fleet based on BN. Before establishing BN, the basic situation of the failure event must be clearly defined as the condition to be studied, usually the potential unsafe condition in the MSAD directive, such as the wing crack event found in this case. As the aircraft's service time accumulates, the wing rib crack grows to fracture, so that the load it bears is transferred to the connected trusses and wing spar webs. The accumulated excess load may cause the skin and wing shell to fail, and ultimately lead to serious consequences such as loss of control of the aircraft.

[0128] When constructing a BN, each node corresponds to a failure propagation-related event one by one, and a directed "edge" is used to represent the causal relationship from the initial cause event to the final unsafe result with a known damage rate. In this case, the BN has a large number of intermediate nodes, and the conditional probability table of each node cannot be obtained in reality. It can be simplified by adding node states and merging intermediate events of the same level. The BN before and after simplification is as follows: Figure 7 As shown, the node information in the figure is as follows Fig. 9As shown in Table 2, in the simplified BN, nodes A to B are binary, "1" represents the failure state, "0" represents the working state; D and E are multi-valued states. The VE method is used to perform causal reasoning on the initial crack failure to solve the posterior probability, that is, the CP value. The specific process is as follows:

[0129] (1) Decompose the BN joint distribution ζ = {P(A), P(BA), P(CB), P(DC), P(ED)};

[0130] (2) Set the evidence variable A = 1 and E as the query variable, and obtain ζ = {P(A = 1), P(B|A = 1), P(C|B), P(D|C), P(E|D)};

[0131] (3) Sort the remaining variables except A and E. Set the elimination order to C, D, B. Eliminate C and we get ζ = {P(A = 1), P(B|A = 1), P(E|D), λ1(B, D)}, where The next variable to be eliminated is D, and we get ζ = {P(A = 1), P(BA = 1), λ2(B, E)}, where the updated factor The last variable to be eliminated is B. Similarly, we can get ζ={P(A=1),λ3(E)}, where

[0132] (4) Multiply and normalize the ζ after elimination to get the posterior probability of the desired variable: Substituting the corresponding conditional probability obtained based on experience, the cumulative probability of all final unsafe results occurring is CP = 0.00364.

[0133] S4. Calculate the total uncorrected fleet (R) based on the events and the result factors. T ), determine its risk level, and output quantitative risk assessment results. The overall uncorrected fleet risk defined by TARAM refers to the statistical number of weighted events expected to occur in the remaining life of the affected fleet if no corrective measures are taken for the fleet when potential unsafe conditions are known. T The corresponding risk threshold for commercial passenger aircraft is 0.02. By comparing the fleet risk value with the threshold to determine whether the fleet is in a dangerous condition, the result can guide safety decisions.

[0134] ND (Not Detected) is the average probability of an undetected defect occurring before an unsafe result or situation occurs, based on historical experience. In this case, ND is 0.2. IR (Injury Ratio) represents the injury rate of the result, i.e., the quantitative form of severity. It is obtained by dividing the number of deaths by the total number of people on the aircraft. It is usually based on historical data of transport aircraft accidents. It can be obtained based on the following data provided by the Transport Aircraft Directorate of the Federal Aviation Administration of the United States: Fig.10Determine from Table 3 in the figure; Since the damage rate of each unsafe result is different, the conditional probability value of each consequence can be calculated according to step S3 and substituted into R T =DA×ND×CP×IR, the final risk result is R T =0.035.

[0135] According to the corresponding risk level, it can be judged that the overall uncorrected risk of the transport fleet in this case is higher than the acceptable risk value of 0.02, indicating that the initial cause of the wing crack event affects the operational safety of the fleet. Therefore, it is necessary to formulate a corresponding risk management plan for the fleet to restore it to an acceptable level.

[0136] In the description of this specification, the description with reference to the terms "one embodiment", "example", "specific example", etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any one or more embodiments or examples in a suitable manner.

[0137] The preferred embodiments of the present invention disclosed above are only used to help illustrate the present invention. The preferred embodiments do not describe all the details in detail, nor do they limit the invention to the specific implementation methods described. Obviously, many modifications and changes can be made according to the content of this specification. This specification selects and specifically describes these embodiments in order to better explain the principles and practical applications of the present invention, so that those skilled in the art can understand and use the present invention well. The present invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for assessing the operational risk of civil aircraft for multiple failure modes, characterized in that: include: Considering the characteristics of multiple wear and failure of civil aircraft structural parts, a failure model of aircraft structural parts based on mixed Weibull distribution is established by using the failure data of fleet operation, and the model parameters are estimated based on the EM algorithm optimized by PSO. Execute Monte Carlo simulation algorithm to simulate the number of defective aircraft with multiple wear and tear failures in the fleet according to the failure model; Construct a Bayesian network to perform failure propagation analysis on the initial cause events that lead to unsafe results, and calculate the probability CP of unsafe results under the condition of the initial cause events; According to the predetermined damage rate IR and non-detection rate ND, the fleet risk level R is calculated. T =DA×ND×CP×IR; The construction of the aircraft structural component failure model includes: Using failure data T = [t i ], construct the likelihood function of complete data: Where N is the total amount of sample data, K is the number of Weibull sub-distributions, and Z = [Z ik ],Z ik is a hidden variable, and when the i-th observation comes from the k-th Weibull subdistribution, then Z ik =1, otherwise Z ik =0; θ is the parameter set of the distribution, expressed as θ=(π1, ..., π k ; θ1, ..., θ k ), π k is the kth sub-distribution weighting factor, θ k =(α k , β k ) represents the scale parameter α of the kth sub-distribution k and shape parameter β k ;f(·|θ k ) is the probability density function of the kth Weibull subdistribution; Execute the EM algorithm by giving the failure data T and the current parameter estimate θ (j) Find the conditional probability P(Z|T,θ) of Z (j) ), obtain the expected value Q of the likelihood function of the conditional probability: Taking the Q function as the objective function, the optimization model for mixed Weibull distribution parameter estimation is established as follows: The EM algorithm also includes: Determine the parameter θ (j) The Q function corresponding to the jth iteration value and Q(θ,θ (j) ) is maximized to obtain the j+1th iteration result θ (j+1) , repeat the above operation until convergence; The EM algorithm based on PSO optimization includes: Correspond the D-dimensional position information of each particle in PSO to the parameter θ; When the algorithm is initialized, the update limits of particle position and velocity are determined according to the initial values ​​of each parameter, and a particle population of popsize is randomly generated; The Q function is used as the fitness function of PSO. The particles are updated by comparing the fitness values ​​of each generation of particles, and the particles with the best fitness values ​​are retained. Update the particles to the maximum number of optimization times and output the global optimal particle position as the maximization result of this round of EM iteration.

2. A method for assessing the operation risk of a civil aircraft for multiple failure modes according to claim 1, characterized in that: Monte Carlo simulation algorithm is executed to simulate the number of defective aircraft with multiple wear and tear failures in the fleet according to the failure model. DA includes: Calculate the cumulative probability of the weighting factor Where q0 = 0, k is the number of Weibull subdistributions; Generate a random number r from a uniform distribution between [0, 1] for each aircraft and determine the probability interval [q k-1 ,q k ], and the aircraft is found to obey the kth Weibull sub-distribution failure mode; Perform inverse transformation simulation to generate aircraft failure time Among them, α k , β k are the scale parameter of the kth sub-distribution and the shape parameter of the kth sub-distribution respectively; Statistics of all aircraft in the fleet during the failure prediction period T m quantity, and T m The average value after multiple simulations is taken as the DA value of the number of defective aircraft.

3. The method for assessing the operation risk of a civil aircraft for multiple failure modes according to claim 1, characterized in that: Constructing a Bayesian network to perform failure propagation analysis on the initial events that lead to unsafe results includes: Analyze the failure propagation chain of aircraft structural parts from initial cause events to unsafe consequences; Each level of intermediate failure corresponds to each node event in the BN; Specific symbols are selected to describe the causal relationship between nodes, and a directed acyclic graph is constructed to intuitively describe the development of unsafe results.

4. The method for assessing the operation risk of a civil aircraft for multiple failure modes according to claim 3 is characterized in that: Calculation of the probability of unsafe outcome under the condition of the initial cause event CP includes: The BN joint distribution is decomposed into the product form of factor functions through variable elimination method; Obtain a factor product formula containing only the query variable by eliminating the corresponding variables in each factor function in an elimination order; The initial event is set as the evidence variable, and the unsafe result is set as the query variable. The distribution of each node is updated through variable elimination and the CP value of the unsafe result under the condition of the initial event occurring is calculated.

Images