Risk optimization allocation method based on covariance matrix adaptive evolution

Abstract

The invention provides a risk optimization distribution method based on self-adaptive evolution of a covariance matrix. The risk optimization distribution method comprises the following steps: firstly, calculating a series of parameters related to positioning precision according to a geometrical shape, observation data and a fault mode of a multi-mode GNSS (Global Navigation Satellite System) satellite of a certain epoch; then determining an initial risk allocation result according to a predetermined integrity risk, a continuity risk and a fault subset number; then taking an initial risk allocation result as an initial value, taking pre-allocated integrity risk and continuity risk as limiting conditions, and taking minimum protection level weighted sum of all subsets as an objective function; an integrity and continuity risk allocation process is abstracted into a mathematical optimization model, adaptive evolution is carried out by adopting a covariance matrix, iterative solution is carried out through processes of sampling, mean value updating, step length control, covariance matrix estimation and the like, and finally an optimal solution of the problem is obtained. Compared with the traditional algorithm, the integrity and continuity risk distributed by the optimization algorithm can significantly reduce the ARAIM protection level.

Description

Technical Field

[0001] The present invention relates to the field of calculation, estimation or counting technology, and in particular to a risk optimization allocation method based on covariance matrix adaptive evolution. Background Art

[0002] As an important auxiliary means for civil aviation approach, GNSS must ensure the quality and performance of its services. Integrity refers to the function of promptly informing users or terminating the signal when the navigation system cannot work according to the predetermined navigation requirements. When GNSS is used for civil aviation approach, integrity plays a more important role than accuracy. The ICAO's standards and recommendations have established strict standards for GNSS navigation performance for terminal approach, which mainly include four performance indicators: accuracy, integrity, continuity and availability. The technology that provides accurate position and descent guidance for aircraft during the landing phase through equipment is called terminal approach technology. Faced with the high safety requirements of civil aviation and the vulnerability of the global navigation satellite system, how to ensure the accuracy, integrity, continuity and availability of GNSS to provide high-quality services for civil aviation has become an urgent problem to be solved.

[0003] In 2006, the Federal Aviation Administration proposed Advanced Receiver Autonomous Integrity Monitoring (ARAIM), which can guide aircraft to an altitude of 200 feet above the ground according to the required navigation performance. Compared with satellite-based augmentation and ground-based augmentation technology, ARAIM does not require a large amount of infrastructure, has low cost and wide coverage, and can provide precision approach services to users in various regions (including polar regions). In addition, because ARAIM uses a multi-hypothesis solution separation algorithm, it can easily and effectively detect a variety of potential satellite failure modes in a multi-mode dual-frequency GNSS environment.

[0004] The traditional ARAIM algorithm evenly distributes continuity risk and integrity risk to the fault subsets of each mode, while ignoring the satellite's own missed alarm rate and prior failure probability, and taking the maximum protection level in the subset as the protection level at the current moment. Although this method is easy to calculate, it cannot guarantee the optimal integrity and continuity risk allocation scheme, resulting in a large protection level and reduced availability. The common methods to improve availability are to increase the number of visible satellites and to improve the user algorithm. Among them, improving the risk optimization allocation method to improve availability is a major focus of current research. The later binary search algorithm did not consider the optimal allocation of continuity risk, and some coefficient optimization methods were prone to singular values ​​when solving, resulting in the inability to solve the protection level, affecting the availability of the algorithm. Summary of the invention

[0005] The present invention aims to address the deficiencies of the prior art in the GNSS risk allocation process and provides a risk optimization allocation method based on the adaptive evolution of the covariance matrix. The method can take into account the allocation of integrity and continuity risks at the same time, and overcome the shortcomings of conservative risk allocation and the appearance of singular values ​​after allocation by finding the optimal solution through a heuristic algorithm. It can better allocate risks and reduce the protection level.

[0006] The present invention provides a risk optimization allocation method based on covariance matrix adaptive evolution, comprising the following steps:

[0007] S1. Obtain a priori failure probability based on GNSS integrity risk support information and failure types to determine the maximum number of satellites or constellations that need to be monitored that fail simultaneously and set target values ​​for GNSS integrity risk and continuity risk;

[0008] S2. Obtain tropospheric noise, multipath noise and observation error noise according to the altitude angle and frequency information, and obtain a weight matrix for evaluating accuracy and continuity, a covariance matrix, a weight matrix for integrity and a covariance matrix corresponding to the weight matrix for integrity according to the tropospheric noise, multipath noise and observation error noise;

[0009] S3, determining the initial risk allocation result through a separation test according to the predetermined integrity risk, continuity risk and number of fault subsets, if the separation test fails, then performing fault detection and elimination and entering step S4; if the separation test passes, then entering step S4;

[0010] S4. Calculate the vertical protection level, take the initial risk allocation result as the initial value, take the pre-allocated integrity risk and continuity risk as the constraint conditions, take the minimum weighted sum of all subset protection levels as the objective function, and convert the risk optimization allocation solution VPL into the solution of the missed alarm rate and false alarm rate of the fault mode;

[0011] Optimization target minY k for:

[0012]

[0013] Among them, st is the restriction condition, VPL k is the vertical protection level of fault subset k, σ v,k is the standard deviation of the fault subset k in the vertical direction, σ ss,k is the standard deviation of the difference between the fault-free subset and the faulty subset k, N set is the total number of fault subsets, K fa,k The false alarm rate P fa,k is the quantile of the standard Gaussian distribution corresponding to the significant level, K md,k The false alarm rate P HMI,kis the quantile of the standard Gaussian distribution corresponding to the significant level, S v,k,r is the weighted least squares matrix of the vertical direction of the rth satellite in subset k, b nom,r is the maximum pseudorange deviation of a single star r under normal conditions, N sat VPL is the number of all visible satellites including both faulty and healthy satellites. k is the vertical protection level of subset k;

[0014] S5. HMI,k and P fa,k As an individual in the population, the optimization target is used as the fitness function, and the covariance matrix is ​​used to adaptively evolve to solve the optimization target minY k The model calculates the final vertical protection level, and a risk optimization allocation method based on the adaptive evolution of the covariance matrix is ​​used to complete it.

[0015] The risk optimization allocation method based on adaptive evolution of covariance matrix described in the present invention, as a preferred embodiment, step S5 comprises the following steps:

[0016] S51, set the covariance matrix adaptive evolutionary algorithm parameters, set the maximum number of iterations, set the number of offspring individuals λ, the number of parent individuals μ (μ<λ), and the reorganization weight w i＝1,...,μ , set the constant value required for adaptation;

[0017] S52, according to the number of subsets N set , in [0,P HMI ] and [0,P fa ] Each generates λ groups of random numbers to form the initial population Z:

[0018]

[0019] in,

[0020] S53, population sampling, weighted recombination of the current optimal subgroup to obtain a new distribution mean:

[0021] and

[0022] in, is the i-th individual in the t+1 generation population, is the average moving step length of the t+1 generation population, m (t) is the mean of the population distribution, σ (t) is the search step length, C (t) is the covariance matrix of the population distribution, N(0, C (t) ) has a mean of 0 and a covariance matrix of C (t) Multivariate Gaussian distribution of ;

[0023] S54, evaluation and selection, calculate the fitness minY of the λ offspring of the contemporary population k , and then sort them according to their fitness, and select the smallest μ individuals Form the best subgroup as the parent sample of the next generation;

[0024] S55, update the mean, and update the strategy parameters by using the information of the optimal subgroup consisting of μ individuals generated by the parent generation: mean m (t) , step size σ (t) , the covariance matrix C (t) , retaining the parent's information:

[0025]

[0026] in, represents the individual ranked i in the t+1 generation population, and the ranking is based on the objective function f = minY for each individual k Determine from small to large;

[0027] S56. Adaptive covariance matrix adjustment: constructing evolutionary paths using exponential smoothing

[0028]

[0029] in, Choose the quality for variance efficiency, is the evolutionary path of the tth generation;

[0030] Combining the rank 1 of the covariance matrix and the rank μ of the covariance matrix, update the covariance matrix:

[0031]

[0032] in, is the covariance evolution path of the t+1th generation, is the transposed matrix of the covariance evolutionary path of the t+1th generation;

[0033] S57, global step length control, Compare with the normalized expected step size E|N(0,I)|:

[0034]

[0035] Among them, exp() is the natural exponential function, under the condition of stationarity That is, the search path can be regarded as a random vector of n-dimensional standard normal distribution, I is the unit matrix, N(0, I) is a multivariate Gaussian distribution with mean 0 and covariance as the unit covariance matrix;

[0036] Compute the covariance evolution path:

[0037]

[0038] S58, until the maximum number of iterations is reached. If so, stop, select the optimal risk allocation scheme according to the optimization goal in the population of the last iteration, and calculate the final vertical protection level; otherwise, return to step S53.

[0039] The risk optimization allocation method based on adaptive evolution of covariance matrix described in the present invention is preferably, in step S51, the constant value includes: the cumulative learning rate α controlled by the step size σ , the learning rate α of the covariance matrix rank 1 c1 , the learning rate α of the covariance matrix rank μ cμ , the evolutionary path decay rate α cp , sampling factor d σ .

[0040] The risk optimization allocation method based on adaptive evolution of covariance matrix described in the present invention is a preferred method. In step S1, the GNSS integrity risk support information includes: the single satellite priori failure probability P within a specified time sat，r , the probability of constellation failure within a specified time P const,j , Space signal accuracy / user ranging accuracyσ URA , spatial signal error / user ranging errorσ URE and the maximum pseudorange deviation b under normal conditions used to assess integrity nom ;

[0041] The prior failure probability P of all satellites sat，r When the satellite prior failure probability is equal A priori probability of failure of constellation for:

[0042]

[0043]

[0044] Among them, N and M are the number of satellites that can be tracked and the number of constellations to which the satellites belong, n s and n c are the number of satellites and constellations assumed to have simultaneous failures, respectively; The number of combinations of the number of trackable satellites to be taken out from the satellites that fail simultaneously; The number of combinations of the number of constellations to which trackable satellites belong from the constellations that fail simultaneously.

[0045] The risk optimization allocation method based on adaptive evolution of covariance matrix described in the present invention is a preferred mode. When the priori failure probability of the failure model is determined, the maximum number of simultaneous failures of satellites or constellations that need to be monitored is also determined.

[0046] The risk optimization allocation method based on adaptive evolution of covariance matrix described in the present invention, as a preferred embodiment, step S3 includes the following steps:

[0047] S31. Constructing a solution separation test quantity: Calculate the least squares location solution for each subset of each fault mode, and construct a test value for solution separation test by subtracting the fault subset solution from the fault-free subset solution:

[0048]

[0049] in, is the location solution of fault subset k in the q direction, is the location solution in the direction of the fault-free subset q, where q is 1, 2 or 3;

[0050] S32, solution separation test threshold calculation: For each fault mode, there are three test thresholds corresponding to the east direction, north direction and zenith direction. The test threshold in the q direction of the fault subset k is T k,q :

[0051] T v,k =K fa,q σ ss,k ;

[0052] Among them, σ ss,k K is the standard deviation of the difference between the fault-free subset and the faulty subset k. fa,q P is the continuous risk fa The significance level corresponds to the standard Gaussian distribution 1-P fa / 2 quantile;

[0053] S33, solution separation test: For all subsets k and directions q, the condition for the solution separation test to pass is:

[0054]

[0055] All subsets of each failure mode are tested separately;

[0056] S34. If the separation test fails, the process proceeds to step S4 after fault detection and elimination; if the classification test passes, the process proceeds to step S4.

[0057] The risk optimization allocation method based on adaptive evolution of covariance matrix described in the present invention is, as a preferred embodiment, in step S4, the vertical protection level is:

[0058] VPL=max(VPL0,max(VPL k ));

[0059] Among them, VPL0 is the vertical protection level of the fault-free subset, and VPL takes the maximum value of the protection levels corresponding to all subsets.

[0060] The risk optimization allocation method based on adaptive evolution of covariance matrix described in the present invention is a preferred method.

[0061]

[0062] Among them, S 0,v,r is the weighted least squares matrix of the rth satellite in the vertical direction of the fault-free subset, S v,k,r is the weighted least squares matrix of the vertical direction of the rth satellite in subset k, T v,k is the test threshold of fault subset k in the vertical direction, σ v,0 is the standard deviation of the fault-free subset in the vertical direction, K md,0 Integrity risk P HMI,0 is the standard Gaussian distribution corresponding to the significant level 1-P HMI,0 / 2 quantile.

[0063] The present invention first calculates a series of parameters related to positioning accuracy based on the geometry, observation data and failure mode of a multi-mode GNSS satellite in a certain epoch; then determines the initial risk allocation result based on the predetermined integrity risk, continuity risk and number of failure subsets; then uses the initial risk allocation result as the initial value, the pre-allocated integrity risk and continuity risk as the constraint condition, and the weighted minimum of all subset protection levels as the objective function; abstracts the integrity and continuity risk allocation process into a mathematical optimization model, adopts the adaptive evolution of the covariance matrix, and finally obtains the optimal solution to the problem through sampling, updating the mean, controlling the step size, and estimating the covariance matrix. The integrity and continuity risks allocated by this optimization algorithm can significantly reduce the ARAIM protection level compared with the traditional algorithm.

[0064] The present invention has the following advantages:

[0065] (1) Different from the traditional optimization algorithm, the present invention uses the adaptive evolution of the covariance matrix and can optimize the objective function without using gradients. It is mainly used to solve difficult nonlinear, non-convex, and black-box optimization problems in continuous domains. It has strong global optimization capabilities, wide applicability, flexible use, and high operating efficiency. It can overcome the shortcomings of the traditional method of conservative risk allocation, which leads to a reduction in the protection level and thus the availability of the ARAIM algorithm.

[0066] (2) The decision value of the traditional method may be close to infinitesimal during the optimization iteration process, resulting in singular values ​​when taking derivatives of infinitesimal variables. The simulated annealing algorithm based on single-objective optimization can avoid this situation and improve the availability of the ARAIM algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 Flowchart of a risk optimization allocation method based on adaptive evolution of covariance matrix. DETAILED DESCRIPTION

[0068] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

[0069] Example 1

[0070] like Figure 1 As shown, a risk optimization allocation method based on adaptive evolution of covariance matrix includes the following steps:

[0071] According to the integrity risk and continuity risk pre-defined in LPV-200, the total integrity risk and continuity risk are pre-allocated to the integrity risk P in the vertical direction. HMI Set to 9.8×10 -8 , continuous risk P fa Set to 3.9×10 -6 ;

[0072] Determine the failure mode based on the pre-set integrity support information and relevant requirements of LPV-200, and calculate the prior failure probability of the corresponding failure mode;

[0073] Integrity support information is a set of information that is pre-set based on the commitment of the constellation manager and the historical performance of the constellation operation. It is used to characterize and express the service quality of the satellite. It includes the following contents.

[0074] Integrity support information content

[0075]

[0076]

[0077] Assuming that all satellites have the same probability of failure, the prior failure probability P of failure mode r is ap,r The calculation is as follows:

[0078]

[0079]

[0080] Among them, N and M are the number of satellites that can be tracked and the number of constellations to which the satellites belong, P sat and P const is the single-satellite failure probability and constellation failure probability in the integrity support information, n s and n c are the number of satellites and constellations assumed to have simultaneous failures, respectively. is the number of combinations of taking a objects out of b objects, Where ! is the factorial. For a positive integer n, its factorial n! is the product of all positive integers from 1 to n, that is, n! = 1×2×3×...×n. When the priori failure probability of the fault model is determined, the maximum number of satellites or constellations that need to be monitored and have simultaneous failures is also determined.

[0081] The simulation experiment uses precision products provided by the IGS Analysis Center of Wuhan University. The ground simulation area is the entire world, with a grid of 10°×10°, a simulation day, a calculation interval of 15 minutes, and a cutoff elevation angle of 10°. The satellite position and ground point position obtained from the precise ephemeris are linearized to obtain the error equation of pseudo-range single-point positioning in the case of a fault-free satellite, which is used to solve the positioning solution calculation of the subset of the separation test process.

[0082] Then go through the following four steps:

[0083] (1) Calculate the integrity, accuracy and continuity weight matrix

[0084] Before calculating the weight matrix, the noise of the troposphere, multipath and observation errors should be calculated based on the altitude angle, frequency information and model.

[0085] According to the broadcast ephemeris and observation files, the satellites with poor observation geometry are eliminated according to the altitude cutoff angle of 10°, and N visible satellites are obtained. The linearized pseudo-range observation equation of all visible satellites without faults and its least square form can be expressed as:

[0086] y = Gx + Δ;

[0087] Where y is the pseudorange residual vector, x is the parameter vector to be estimated, including the receiver position parameters and clock error parameters, G is the coordinate coefficient matrix, and Δ is the pseudorange error vector.

[0088] Before calculating the weight matrix, the tropospheric noise σ should be calculated based on the altitude angle and frequency information. n,tropo , multipath noise σ MP .

[0089]

[0090]

[0091]

[0092]

[0093] Where θ is the altitude angle, f1 and f2 are the two frequencies used in the deionosphere combination, and σ noise is the observation noise, σ n,user For the comprehensive impact on users.

[0094] Using the above model, we can get the error noise, and according to the noise, we can calculate the weight matrix W for evaluating accuracy and continuity. URE and the corresponding covariance matrix C acc , the integrity weight matrix W URA and the corresponding covariance matrix C int .

[0095]

[0096]

[0097] Among them, URA is the user ranging accuracy parameter.

[0098] The calculated weight matrix and covariance matrix will be used to locate the solution of the calculation subset in step (2).

[0099] (2) Constructing the solution separation test volume

[0100] The solution separation test is to calculate the least squares location solution of each subset of each fault mode, and then construct the test value Δx by subtracting the faulty subset solution from the non-faulty subset solution. k , used for solution separation test.

[0101]

[0102] In the formula, represents the location solution in the q direction of the fault subset k, represents the location solution in the direction of the fault-free subset q.

[0103] (3) Solution separation test threshold calculation

[0104] For each fault mode, there are three test thresholds corresponding to the east direction, north direction and zenith direction. For the fault subset k, the test threshold in the direction of q can be expressed as T k,q .

[0105] T v,k =K fa,q σ ss,k ;

[0106] σ ss,kK is the standard deviation of the difference between the fault-free subset and the faulty subset k. fa,q P is the continuous risk fa is the standard Gaussian distribution corresponding to the significant level quantile of .

[0107] (4) Solution separation test

[0108] For all subsets k and directions q, the solution separation test passes under the following conditions:

[0109]

[0110] All subsets of each fault mode are subjected to separation test. Only when the separation test is passed can the corresponding protection level be calculated. If the condition is not met, fault detection and elimination should be carried out.

[0111] The vertical protection level can be expressed as follows:

[0112] VPL=max(VPL0,max(VPL k ));

[0113] In the above formula, VPL0 represents the vertical protection level of the fault-free subset, VPL k It represents the vertical protection level of fault subset k. Finally, VPL takes the maximum value of the protection levels corresponding to all subsets.

[0114]

[0115]

[0116] Where N sat is the number of all visible satellites; S 0,v,r is the weighted least squares matrix of the rth satellite in the vertical direction of the fault-free subset; S v,k,r is the weighted least squares matrix of the vertical direction of the rth satellite in subset k; T v,k is the test threshold of fault subset k in the vertical direction, where σ v,0 is the standard deviation of the fault-free subset in the vertical direction, σ v,k is the standard deviation of the fault subset k in the vertical direction; K md,0 and K md,k The integrity risk P HMI,0 and P HMI,k is the standard Gaussian distribution corresponding to the significant level and Therefore, the calculation of VPL is transformed into the optimization of the distribution of integrity risk P HMI and continuity risk P fa, solve the allocation scheme with the largest subset and the smallest VPL among all allocation schemes, and use the VPL calculated by this scheme as the final result.

[0117] The optimal solution of the problem is obtained by defining the optimization problem and setting constraints. Based on the protection level calculation model and constraints, the optimization goal can be abstracted as follows:

[0118] min(max(VPL k (P HMI.k ,P fa,k ))

[0119]

[0120] Where P HMI,k is the integrity risk assigned to the fault subset k; P fa,k is the continuity risk assigned to the fault subset k; P fa For the continuous risk pre-assigned to the vertical direction; P HMI The integrity risk pre-assigned to the vertical direction. k (P HMI.k ,P fa,k ) is the vertical protection level of subset k, max() is the maximum value function, and min() is the minimum value function. st is the constraint condition. The optimization goal is to solve the risk allocation result of the subset with the largest vertical protection level in all allocation schemes under the constraint that the sum of the integrity and continuity risks of all subsets is less than the pre-allocated integrity and continuity risks in the vertical direction through the risk allocation algorithm. In general, the VPL corresponding to each fault subset is equal, so the optimization goal can be rewritten as follows:

[0121] VPL = min(VPL k )k={0,···,N set};

[0122] According to multi-objective optimization theory, multi-objective optimization needs to be transformed into a single-objective optimization problem. k Weighted sum as optimization target minY k

[0123]

[0124] P HMI,k and P fa,k As individuals in the population, the optimization target is used as the fitness function, and the covariance matrix is ​​used to adaptively evolve to solve the optimization model. The specific algorithm steps are as follows:

[0125] Step 1: Set the covariance matrix adaptive evolutionary algorithm parameters, set the maximum number of iterations, set the number of offspring individuals λ, the number of parent individuals μ (μ<λ), and the reorganization weight w i＝1,...,μ , set the constant values ​​required for adaptation, including the cumulative learning rate α controlled by the step size σ , the learning rate α of the covariance matrix rank 1 c1 , the learning rate α of the covariance matrix rank μ cμ , the evolutionary path decay rate α cp , sampling factor d σ ;

[0126] Step 2: According to the number of subsets N set , in [0,P HMI ] and [0,P fa ] Each generates λ groups of random numbers to form the initial population Z:

[0127]

[0128] Among them are

[0129] Step 3: Population sampling, weighted recombination of the current optimal subgroup to obtain a new distribution mean.

[0130]

[0131] is the i-th individual in the t+1th generation population, is the average moving step length of the t+1 generation population, m (t) is the mean of the population distribution, σ (t) is the standard deviation of the population distribution (search step length), C (t) is the covariance matrix of the population distribution;

[0132] Step 4: Evaluation and selection, calculate the fitness minY of the λ offspring of the contemporary population k , and then sort them according to their fitness, and select the smallest μ individuals Form the optimal subgroup as the parent sample of the next generation.

[0133] Step 5: Update the mean. The algorithm updates the strategy parameters by using the information of the optimal subgroup consisting of μ individuals generated by the parent generation: mean m (t) , step size σ (t) , the covariance matrix C (t) Keep the parent's information.

[0134]

[0135] represents the individual ranked i in the t+1 generation population, and the ranking is based on the objective function f = minY for each individual k Determine from small to large.

[0136] Step 6: Adaptive covariance matrix adjustment. In actual calculation, exponential smoothing is used to construct the evolution path. as follows:

[0137]

[0138] in, Choose quality for variance efficiency.

[0139] Then, combining the rank 1 of the covariance matrix and the rank μ of the covariance matrix, update the covariance matrix:

[0140]

[0141] Step7: Global step length control, Compared with its normalized expected step size E|N(0,I)|, we have:

[0142]

[0143] exp() represents the natural exponential function, under the condition of stationarity That is, the search path can be regarded as a random vector with an n-dimensional standard normal distribution.

[0144] Then, calculate the covariance evolution path:

[0145]

[0146] Step 8: Until the maximum number of iterations is reached, if it is reached, stop, select the optimal risk allocation plan according to the optimization goal in the population of the last iteration, and calculate the final vertical protection level, otherwise return to Step 3.

[0147] The global simulated vertical protection level calculated after the optimal allocation of GNSS integrity and continuity risks based on the covariance matrix adaptive evolution is compared with two traditional algorithms (from left to right, they are the risk average allocation algorithm, the binary search algorithm, and the covariance matrix adaptive evolution algorithm). The protection level obtained by the covariance matrix adaptive evolution algorithm is reduced to varying degrees in different regions of the world, which is 33.9% and 8.9% lower than that of the risk average allocation algorithm and the binary search algorithm, respectively, which also verifies the effectiveness of the method.

[0148] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical scheme and inventive concept of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. A risk optimization allocation method based on adaptive evolution of covariance matrix, characterized by: The following steps are involved: S1. Obtain a priori failure probability based on GNSS integrity risk support information and failure types to determine the maximum number of satellites or constellations that need to be monitored that fail simultaneously and set target values ​​for GNSS integrity risk and continuity risk; S2. Obtain tropospheric noise, multipath noise and observation error noise according to the altitude angle and frequency information, and obtain a weight matrix for evaluating accuracy and continuity, a covariance matrix, a weight matrix for integrity and a covariance matrix corresponding to the weight matrix for integrity according to the tropospheric noise, the multipath noise and the observation error noise; S3, determining the initial risk allocation result through a separation test according to the predetermined integrity risk, continuity risk and number of fault subsets, if the separation test fails, then performing fault detection and elimination and entering step S4; if the separation test passes, then entering step S4; S4. Calculate the vertical protection level, take the initial risk allocation result as the initial value, take the pre-allocated integrity risk and continuity risk as the constraint conditions, take the minimum weighted sum of all subset protection levels as the objective function, and convert the risk optimization allocation solution VPL into the solution of the missed alarm rate and false alarm rate of the fault mode; Optimization target minY k for: Among them, st is the restriction condition, VPL k is the vertical protection level of fault subset k, σ v,k is the standard deviation of the fault subset k in the vertical direction, σ ss,k is the standard deviation of the difference between the fault-free subset and the faulty subset k, N set is the total number of fault subsets, K fa,k The false alarm rate P fa,k is the quantile of the standard Gaussian distribution corresponding to the significant level, K md,k The false alarm rate P HMI,k is the quantile of the standard Gaussian distribution corresponding to the significant level, S v,k,r is the weighted least squares matrix of the vertical direction of the rth satellite in subset k, b nom,r is the maximum pseudorange deviation of a single star r under normal conditions, N sat is the number of all visible satellites including both faulty and non-faulty satellites; S5. HMI,k and P fa,k As an individual in the population, the optimization target is used as the fitness function, and the covariance matrix is ​​used to adaptively evolve to solve the optimization target minY k The model calculates the final vertical protection level, which is completed by a risk optimization allocation method based on the adaptive evolution of the covariance matrix.

2. A risk optimization allocation method based on covariance matrix adaptive evolution according to claim 1, characterized in that: Step S5 includes the following steps: S51, set the covariance matrix adaptive evolutionary algorithm parameters, set the maximum number of iterations, set the number of offspring individuals λ, the number of parent individuals μ (μ<λ), and the reorganization weight w i＝1,...,μ , set the constant value required for adaptation; S52, according to the number of subsets N set , in [0,P HMI ] and [0,P fa ] Each generates λ groups of random numbers to form the initial population Z: in, S53, population sampling, weighted recombination of the current optimal subgroup to obtain a new distribution mean: and in, is the i-th individual in the t+1 generation population, is the average moving step length of the t+1 generation population, m (t) is the mean of the population distribution, σ (t) is the search step length, C (t) is the covariance matrix of the population distribution, N(0, C (t) ) has a mean of 0 and a covariance matrix of C (t) Multivariate Gaussian distribution of ; S54, evaluation and selection, calculate the fitness minY of the λ offspring of the contemporary population k , and then sort them according to their fitness, and select the smallest μ individuals Form the best subgroup as the parent sample of the next generation; S55, update the mean, and update the strategy parameters by using the information of the optimal subgroup consisting of μ individuals generated by the parent generation: mean m (t) , step size σ (t) , the covariance matrix C (t) , retaining the parent's information: in, represents the individual ranked i in the t+1 generation population, and the ranking is based on the objective function f = minY for each individual k Determine from small to large; S56. Adaptive covariance matrix adjustment: constructing evolutionary paths using exponential smoothing in, Choose the quality for variance efficiency, is the evolutionary path of the tth generation; Combining the rank 1 of the covariance matrix and the rank μ of the covariance matrix, update the covariance matrix: in, is the covariance evolution path of the t+1th generation, is the transposed matrix of the covariance evolutionary path of the t+1th generation; S57, global step length control, Compare with the normalized expected step size E|N(0,I)|: Among them, exp() is the natural exponential function, under the condition of stationarity That is, the search path can be regarded as a random vector of n-dimensional standard normal distribution, I is the unit matrix, N(0, I) is a multivariate Gaussian distribution with mean 0 and covariance as the unit covariance matrix; Compute the covariance evolution path: S58, until the maximum number of iterations is reached. If so, stop, select the optimal risk allocation scheme according to the optimization goal in the population of the last iteration, and calculate the final vertical protection level; otherwise, return to step S53.

3. A risk optimization allocation method based on covariance matrix adaptive evolution according to claim 2, characterized in that: In step S51, the constant value includes: the cumulative learning rate α of the step size control σ , the learning rate α of the covariance matrix rank 1 c1 , the learning rate α of the covariance matrix rank μ cμ , evolutionary path decay rate α cp , sampling factor d σ .

4. A risk optimization allocation method based on covariance matrix adaptive evolution according to claim 1, characterized in that: In step S1, the GNSS integrity risk support information includes: the priori failure probability P of a single satellite within a specified time sat，r , the probability of constellation failure within a specified time P const,j , Space signal accuracy / user ranging accuracyσ URA , spatial signal error / user ranging errorσ URE and the maximum pseudorange deviation b under normal conditions used to assess integrity nom ; The prior failure probability P of all satellites sat，r When the satellite prior failure probability is equal A priori probability of failure of constellation for: Among them, N and M are the number of satellites that can be tracked and the number of constellations to which the satellites belong, n s and n c are the number of satellites and constellations assumed to have simultaneous failures, respectively; The number of combinations of the number of trackable satellites to be taken out from the satellites that fail simultaneously; The number of combinations of the number of constellations to which trackable satellites belong from the constellations that fail simultaneously.

5. A risk optimization allocation method based on covariance matrix adaptive evolution according to claim 4, characterized in that: When the priori failure probability of the fault model is determined, the maximum number of satellites or constellations that need to be monitored and have simultaneous failures is also determined.

6. A risk optimization allocation method based on covariance matrix adaptive evolution according to claim 1, characterized in that: Step S3 includes the following steps: S31. Constructing a solution separation test quantity: Calculate the least squares location solution for each subset of each fault mode, and construct a test value for solution separation test by subtracting the fault subset solution from the fault-free subset solution: in, is the location solution of fault subset k in the q direction, is the location solution in the direction of the fault-free subset q, where q is 1, 2 or 3; S32, solution separation test threshold calculation: For each fault mode, there are three test thresholds corresponding to the east direction, north direction and zenith direction. The test threshold in the q direction of the fault subset k is T k,q : T v,k =K fa,q s ss,k ; Among them, σ ss,k K is the standard deviation of the difference between the fault-free subset and the faulty subset k. fa,q P is the continuous risk fa The significance level corresponds to the standard Gaussian distribution 1-P fa / 2 quantile; S33, solution separation test: For all subsets k and directions q, the condition for the solution separation test to pass is: All subsets of each failure mode are tested for separation; S34. If the separation test fails, the process proceeds to step S4 after fault detection and elimination; if the classification test passes, the process proceeds to step S4.

7. A risk optimization allocation method based on covariance matrix adaptive evolution according to claim 1, characterized in that: In step S4, the vertical protection level is: VPL=max(VPL0,max(VPL k )); Among them, VPL0 is the vertical protection level of the fault-free subset, and VPL takes the maximum protection level corresponding to all subsets.

8. A risk optimization allocation method based on covariance matrix adaptive evolution according to claim 7, characterized in that: Among them, S 0,v,r is the weighted least squares matrix of the rth satellite in the vertical direction of the fault-free subset, S v,k,r is the weighted least squares matrix of the vertical direction of the rth satellite in subset k, T v,k is the test threshold of fault subset k in the vertical direction, σ v,0 is the standard deviation of the fault-free subset in the vertical direction, K md,0 Integrity risk P HMI,0 is the standard Gaussian distribution corresponding to the significant level 1-P HMI,0 / 2 quantile.

Images