Patent Publication Number: US-2023137147-A1

Title: Constellation configuration optimization method of leo satellite augmentation system for araim application

Description:
CROSS-REFERENCE TO RELATED APPLICATION 
     This application is a continuation application of International Patent Application No. PCT/CN2021/090796, filed on Apr. 29, 2021, which itself claims priority to and benefit of Chinese Patent Application No. 202110306891.3 filed on Mar. 23, 2021 in the State Intellectual Property Office of P. R. China. The disclosure of each of the above applications is incorporated herein by reference in its entirety. 
    
    
     TECHNICAL FIELD 
     The present invention relates to the technical field of satellite navigation, in particular to a constellation configuration optimization method of a low earth orbit (LEO) satellite augmentation system for an ARAIM application. 
     BACKGROUND ART 
     A Global Navigation Satellite System (GNSS) has developed rapidly and has been widely used in various fields. In the future, mass users all over the world have put forward a series of requirements for high precision, fast convergence, high integrity, high security and high availability, so high-security real-time precise positioning in complex environment has become a higher goal of GNSS at present. In order to solve the above problems, all countries are actively building their own navigation enhancement systems, while the traditional navigation enhancement systems are all working systems of ground station auxiliary information enhancement. However, due to the limitation of our country&#39;s territory and territory, it is difficult to realize global centimeter-level positioning and fast convergence service by setting up global ground stations. Therefore, with the increasing demand of station density and information rate, the existing augmentation system can&#39;t support the next generation Beidou to realize centimeter-level real-time positioning service in challenging environment. Therefore, it is necessary to design the architecture and system of a new enhanced system to achieve three goals: 1) to reduce the dependence on overseas stations; 2) to expand enhanced services from regional enhancement to global enhancement; 3) to greatly improve the usability and security of precision positioning services. 
     With the development of communication services, at the end of the 20th century, LEO satellite constellations for mobile communication appeared, and the typical representatives are Iridium and Globalstar constellations in the United States. In January 2019, the second generation iridium system completed networking, completely replacing the first generation iridium system. In addition to the original communication service, the new generation Iridium satellite system is equipped with navigation enhanced payload, which can provide users with other services besides communication, such as positioning and tracking, broadcast automatic correlation monitoring and so on. Since 2015, many internationally renowned companies, such as OneWeb, SpaceX, Boeing of the United States, Samsung of South Korea, Aerospace Science and Technology Group of China and Aerospace Science and Technology Group, have successively announced the launch and deployment of their commercial LEO constellations, all in order to provide seamless and stable broadband Internet communication services to the world. The first star of China&#39;s “Hongyan” constellation was successfully launched on Dec. 29, 2018 and entered the scheduled orbit. The whole constellation is planned to be deployed around 2024. Hongyan Mobile Communication Constellation not only meets the needs of multi-field monitoring data information transmission, but also has the mobile broadcasting function, is equipped with onboard GNSS receiver, and is equipped with navigation enhancement function. 
     Because of the rapid development of LEO satellite market and the low cost of carrying navigation enhancement payload, it can be used as a space-based monitoring platform to realize a new navigation enhancement system of “space-based monitoring+signal enhancement”. The opening of Beidou-3 marks the official completion of the meter-level navigation and positioning system comparable to GPS in China. For the future development of Beidou, Ran Chengqi, deputy chief designer of Beidou satellite navigation system engineering, said that it is expected to be referred to as the space-based low-orbit constellation system for short before 2025, and provide centimeter-level positioning services to the whole world. 
     At present, the research on LEO navigation enhancement mainly focuses on real-time precise orbit determination and remote sensing monitoring. There is a lack of integrity monitoring of LEO satellite navigation enhancement system, and once LEO satellite navigation system is officially used for navigation enhancement service, it will affect the availability of ARAIM airborne receiver. 
     Therefore, in order to overcome the problems in the prior art, a constellation configuration optimization method of a LEO satellite augmentation system for an ARAIM application is needed. 
     SUMMARY 
     An object of the present invention is to provide a constellation configuration optimization method of a LEO satellite augmentation system for an ARAIM application, which includes the following steps: 
     1, traversing vertical protection levels after all subset solutions and fault modes under the condition that integrity risk and continuity risk are equally distributed, and determining the constraint conditions of LEO satellite constellation configuration parameters; 
     2, determining objective functions of LEO satellite constellation configuration parameters x 1 , x 2 , x 3 , x 4 , eliminating calculated values of abnormal vertical protection levels, and screening initial populations of the parameters x 1 , x 2 , x 3 , x 4 ; 
     3, calculating fitness of the objective functions of the LEO satellite constellation configuration parameters x 1 , x 2 , x 3 , x 4 ; 
     4, after screening out an initial population of the parameters x 1 , x 2 , x 3 , x 4 , starting from a second generation population, merging a parent population and an offspring population to form a new offspring population, and performing fast non-dominated sorting, calculating crowding degrees of individuals in each non-dominated layer, randomly paring the individuals, and performing a crossover operation of a genetic algorithm between the paired two individuals; 
     5, after the parent population and the offspring population are merged to form the new offspring population, implementing optimal preservation strategy and local optimal selection for the new offspring population, and selecting a maximum value of the objective functions as an optimal offspring, and repeating step  4  until a genetic algebra is less than a maximum genetic algebra. 
     Preferably, the vertical protection level in step 1 is expressed as: 
     VPL q =max((VPL 0 ) q ,max ((VPL i ) q ), where VPL 0  is a vertical protection level under a fault-free condition, (VPL 0 ) q =K MD,q   0 ·σ 0,q +Σ n=1   Nsat |S q,n   0 |·b int,n , K MD,q   0  is an integrity and continuity risk value of the fault-free mode under a fully visible satellite subset, S q,n   0 , is a projection matrix, and b int,n  is a maximum nominal deviation of a n th  satellite; 
     VPL, is a vertical protection level corresponding to a measurement deviation of an i th  fault mode where a maximum deviation is not exceeded, (VPL i ) q =K MD,q   i ·σ i,q +Σ n=1   Nsat |S q,n   i |·b int,n +D i,q , D i,q  is a detection threshold; 
     wherein , σ 0,q =√{square root over ((GW INT G T ) q,q   −1 )}, σ i,q =√{square root over ((G T M i W INT G T ) q,q   −1 )}, G is a geometric matrix in a pseudo-range observation equation, W INT  is a parameter of a fixed error hypothesis model related to integrity, M t  is an identity matrix of a size N sat *N sat , N sat  is the number of visible satellites, and q is a q th  sample point; 
     the geometric matrix G in the pseudo-range observation equation contains the LEO satellite constellation configuration parameters x 1 , x 2 , x 3 , x 4 ; 
     the constraint conditions of the LEO satellite constellation configuration parameters are expressed as: 
     
       
         
           
             
               
                 
                   [ 
                   
                     
                       
                         
                           N 
                           0 
                         
                       
                       
                         0 
                       
                     
                     
                       
                         
                           N 
                           C 
                         
                       
                       
                         
                           N 
                           SO 
                         
                       
                     
                   
                   ] 
                 
                 [ 
                 
                   
                     
                       
                         Ω 
                         ij 
                       
                     
                     
                       
                         Ω 
                         11 
                       
                     
                   
                   
                     
                       
                         M 
                         ij 
                       
                     
                     
                       
                         M 
                         11 
                       
                     
                   
                 
                 ] 
               
               = 
               
                 
                   2 
                   π 
                 
                 [ 
                 
                   
                     
                       i 
                     
                     
                       
                         - 
                         1 
                       
                     
                   
                   
                     
                       j 
                     
                     
                       
                         - 
                         1 
                       
                     
                   
                 
                 ] 
               
             
             , 
           
         
       
     
     where N 0  is the number of orbital planes of LEO constellation, N SO  is the number of satellites on each orbital plane, N C  is a phase parameter, N C ∈[1, N 0 ]. Ω is a right ascension of an ascending intersection point, M is an average anomaly, where i represents an i th  orbital plane and j represents j th  satellite. 
     Preferably, the objective functions of the LEO satellite constellation configuration parameters x 1 , x 2 , x 3 , x 4  in step 2 are expressed as: 
     
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             1 
                           
                           ) 
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         0 
                       
                       ≤ 
                       
                         x 
                         1 
                       
                       ≤ 
                       
                         π 
                         2 
                       
                     
                   
                 
               
               , 
             
           
         
       
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             2 
                           
                           ) 
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         800 
                       
                       ≤ 
                       
                         x 
                         2 
                       
                       ≤ 
                       2000 
                     
                   
                 
               
               , 
             
           
         
       
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             3 
                           
                           ) 
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         0 
                       
                       ≤ 
                       
                         x 
                         3 
                       
                       ≤ 
                       0.782 
                     
                   
                 
               
               , 
             
           
         
       
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             4 
                           
                           ) 
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         0 
                       
                       ≤ 
                       
                         x 
                         4 
                       
                       ≤ 
                       0.3 
                     
                   
                 
               
               , 
             
           
         
       
     
     wherein, F VPL ( x 1 ) is a function where a vertical protection level VPL is the parameter x 1 , F VPL ( x 2 ) is a function where a vertical protection level VPL is the parameter x 2 , F VPL ( x 3 ) is a function where a vertical protection level VPL is the parameter x 3 , and F VPL ( x 4 ) is a function where a vertical protection level VPL is the parameter x 4 . 
     Preferably, in step 2, the calculated values of abnormal vertical protection levels are eliminated by the following method: 
     initial anomaly detection thresholds are set to T first   + , T first   − , and T first   + =μ all +2σ all , T first   − =μ all −2σ all  is defined, where 
     μ all  is an average value of VPL calculated values of vertical protection levels corresponding to all sample points, and σ all  is a standard deviation of vertical protection and VPL calculated values corresponding to all sample points, 
     the calculated value of the vertical protection level VPL of each sample point is compared with T first   + , T first   − , and if the inequality T first   + ≤VPL≤T first   −  is satisfied, the sample point passes the threshold detection and waits for the initial population screening; if the inequality T first   + ≤VPL≤T first   −  is not satisfied, the sample point fails the threshold detection and is stored in a abnormal data module. 
     Preferably, in step 2, the initial population screening is carried out by the following method: 
     for an orbit inclination parameter x 1 , a sampling interval is set to 0.01, and a total of 158 sample points are generated; a set of data is taken every 9 sample points as a sample of the initial population, and the 15 samples constitute the initial population of the orbit inclination parameter x 1 ; 
     for an orbit height parameter x 2 , a sampling interval is set to 1, and a total of 1200 sample points are generated; a set of data is taken every 19 sample points as a sample of the initial population, and the 60 samples form the initial population of the orbit inclination parameter x 2 ; 
     for an initial value parameter x 3  of ascending intersection right ascension, a sampling interval is set to 0.001, and a total of 79 sample points are generated; a set of data is taken every four sample points as a sample of the initial population, and the 15 samples constitute the initial population of the orbit inclination parameter x 3 ; 
     for an initial value parameter x 4  of mean anomaly, a sampling interval is set to 0.01, and a total of 30 sample points are generated; a set of data is taken every four sample points as a sample of the initial population, and the six samples constitute the initial population of the orbit inclination parameter x 4 . 
     Preferably, in step 3, the fitness of the objective functions of the LEO satellite constellation configuration parameters x 1 , x 2 , x 3 , x 4  is calculated by the following method: 
     a maximum optimization problem function is adopted for a fitness function: 
     
       
         
           
             
               Fitness 
               ⁢ 
                   
               
                 ( 
                 
                   
                     F 
                     VPL 
                   
                   ( 
                   x 
                   ) 
                 
                 ) 
               
             
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           
                             
                               F 
                               VPL 
                             
                             ( 
                             x 
                             ) 
                           
                           + 
                           
                             C 
                             min 
                           
                         
                         , 
                         
                           
                             
                               C 
                               min 
                             
                             + 
                             
                               
                                 F 
                                 VPL 
                               
                               ( 
                               x 
                               ) 
                             
                           
                           &gt; 
                           0 
                         
                       
                     
                   
                   
                     
                       
                         0 
                         , 
                         
                           
                             
                               C 
                               min 
                             
                             + 
                             
                               
                                 F 
                                 VPL 
                               
                               ( 
                               x 
                               ) 
                             
                           
                           ≤ 
                           0 
                         
                       
                     
                   
                 
                 , 
               
             
           
         
       
     
     where C min  is a preset number, and the minimum function value of the objective function F VPL (x) estimated so far is taken as the objective function, F VPL (x) is an objective function, representing a function where the vertical protection level VPL is the parameter x 1 , x 2 , x 3  or x 4 . 
     Preferably, the objective function should also satisfy the following conditions: 
     the value of the objective function is ≤35 m. 
     Preferably, the merging ratio of the parent population and the offspring population to form the new offspring population is: 
     
       
         
           
             
               
                 pro 
                 i 
               
               = 
               
                 
                   
                     ❘ 
                     &#34;\[LeftBracketingBar]&#34; 
                   
                   
                     
                       
                         X 
                         ul 
                       
                       - 
                       
                         X 
                         ll 
                       
                     
                     
                       Δπ 
                       i 
                     
                   
                   
                     ❘ 
                     &#34;\[RightBracketingBar]&#34; 
                   
                 
                 
                   2 
                   ⁢ 
                   
                     ( 
                     
                       
                         N 
                         interval 
                       
                       + 
                       1 
                     
                     ) 
                   
                 
               
             
             , 
             
               i 
               = 
               1 
             
             , 
             2 
             , 
             3 
             , 
             4 
             , 
           
         
       
     
     where X ul  is an upper limit of an optimization parameter range, X ll  is a lower limit of the optimization parameter range, Δτ i  is an sampling interval of each parameter, and N interval  is a sample interval when the initial population is generated. 
     Preferably, in step 5, a local mean A m   n  of the sample is introduced for local optimal selection, a selection threshold T sel  is set, and after performing the optimal preservation strategy selection for the new offspring population, a difference between the objective function and each local mean A m   n  is calculated, 
     if the difference is greater than the selection threshold T sel , the local optimal test is passed, and the maximum value of the objective function is taken as the optimal offspring; if the difference is less than the threshold T sel , the maximum value near the objective function will be searched as the optimal offspring. 
     10. The method according to claim  9 , wherein the local mean A m   n  of the sample is expressed as: 
     
       
         
           
             
               
                 A 
                 m 
                 n 
               
               = 
               
                 
                   
                     
                       ∑ 
                       
                         
                           5 
                           ⁢ 
                           n 
                         
                         + 
                         1 
                       
                       
                         5 
                         ⁢ 
                         
                           ( 
                           
                             n 
                             + 
                             1 
                           
                           ) 
                         
                       
                     
                     
                       x 
                       i 
                     
                   
                   k 
                 
                 ⁢ 
                 
                   ( 
                   
                     
                       n 
                       = 
                       0 
                     
                     , 
                     1 
                     , 
                     2 
                     , 
                     3 
                     , 
                     … 
                        
                     , 
                     
                       m 
                       = 
                       1 
                     
                     , 
                     2 
                     , 
                     3 
                     , 
                     … 
                        
                     , 
                     
                       [ 
                       
                         
                           [ 
                           
                             
                               
                                 X 
                                 ul 
                               
                               - 
                               
                                 X 
                                 ll 
                               
                             
                             
                               Δ 
                               ⁢ 
                               
                                 π 
                                 i 
                               
                             
                           
                           ] 
                         
                         5 
                       
                       ] 
                     
                   
                   ) 
                 
               
             
             , 
           
         
       
     
     where m represents the m th  local mean value, the maximum value is 
     
       
         
           
             
               [ 
               
                 
                   [ 
                   
                     
                       
                         X 
                         ul 
                       
                       - 
                       
                         X 
                         ll 
                       
                     
                     
                       Δ 
                       ⁢ 
                       
                         π 
                         i 
                       
                     
                   
                   ] 
                 
                 5 
               
               ] 
             
             , 
           
         
       
     
     and x i  represents the i th  sample of the corresponding parameter, it is stipulated that the average value of every five samples is taken as the local mean value, X ul  is an upper limit of the optimization parameter range, X ll  is a lower limit of the optimization parameter range, and Δτ i  is a X ul  sampling interval of each parameter. 
     According to the present situation of LEO navigation enhancement, the present invention provides a constellation configuration optimization method suitable for a LEO satellite navigation enhancement system, and a constellation configuration optimization method based on an ARAIM protection level algorithm is proposed by combining the non-dominated sorting genetic algorithm with elite strategy. The ARAIM protection level algorithm optimization is realized from another angle through operations such as abnormal data elimination and parameter sampling, which makes up for the integrity monitoring vacancy of LEO satellite navigation enhancement system and provides reference for future LEO satellite navigation system design and networking. 
     According to the method, the parameters to be optimized are determined by the ARAIM protection level formula, and the specific observation matrix to be optimized is obtained after the risk values are evenly distributed; Define the constraint range of the optimization parameters in the observation matrix, and optimize the objective function within the constraint range; before optimization, abnormal data detection and data sampling are carried out to screen the initial population, and when the parent and offspring populations are merged, a specific proportion is adopted to merge them purposefully; Finally, in order to avoid local optimization, the local mean of the sample is defined, and the optimized value and the local mean are threshold-tested, so as to achieve the goal of global optimization. 
     It should be understood that the foregoing general description and the following detailed description are exemplary illustrations and explanations, and should not be used as limitations on the claimed contents of the present invention. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       Further objects, functions and advantages of the present invention will be clarified by the following description of embodiments of the present invention with reference to the accompanying drawings, in which: 
         FIG.  1    schematically shows the spatial schematic diagram of the LEO enhanced ARAIM system of the present invention. 
         FIG.  2    shows the algorithm flow chart of constellation configuration optimization method of LEO augmentation system for ARAIM application. 
         FIG.  3    shows an initial anomaly verification diagram of the vertical protection level. 
     
    
    
     DETAILED DESCRIPTION 
     With reference to exemplary embodiments, the objects and functions of the present invention and the methods for realizing them will be clarified. However, the present invention is not limited to the following disclosed exemplary embodiments; It can be implemented in different forms. The essence of the description is only to help those skilled in relevant fields to comprehensively understand the specific details of the present invention. 
     Hereinafter, embodiments of the present invention will be described with reference to the accompanying drawings. In the drawings, the same reference numerals represent the same or similar parts or the same or similar steps. 
     Advanced Receiver Autonomous Integrity Monitor (ARAIM) is an important technology in satellite navigation integrity monitoring. If there is any fault in the airborne terminal, it can alert the user in time, without building a large number of ground infrastructure. Its application is convenient and quick, and it can be quickly popularized. 
     In order to solve the problems existing in the prior art, the present invention aims to improve the availability of ARAIM of the Beidou satellite navigation system, and reduce its protection level by introducing LEO satellite navigation system. The present invention selects the optimal LEO constellation configuration for different parameters by using non-dominated sorting genetic algorithm with elite strategy, so as to improve the availability of an ARAIM airborne receiver. 
     In order to explain the present invention more clearly, firstly, the LEO satellite system is briefly described. As shown in  FIG.  1   , the LEO enhanced ARAIM system of the present invention is a spatial schematic diagram. The LEO satellite system is equipped with an onboard GNSS receiver, which generates a signal with the same time-frequency as the Beidou satellite navigation system, and broadcasts it to the ground through the onboard navigation load. After receiving the communication navigation signal, the ground monitoring station transmits it to the main control station, and the main control station determines the orbit information of the LEO satellite after calculation, and uploads it to the LEO satellite through the uplink of the injection station. At this time, the LEO satellite can download the ephemeris containing the orbit information to the user terminal through the downlink, and the user terminal (the airborne terminal) can optimize the constellation configuration of the LEO satellite augmentation system by using this information to solve the problem with an algorithm. 
       FIG.  2    shows the algorithm flow chart of the constellation configuration optimization method of a LEO augmentation system for an ARAIM application of the present invention. According to the embodiment of the present invention, a constellation configuration optimization method of LEO satellite augmentation system for ARAIM application includes the following steps: 
     1. Traversing vertical protection levels after all subset solutions and fault modes under the condition that integrity risk and continuity risk are equally distributed, and determining the constraint conditions of LEO satellite constellation configuration parameters; 
     After the ephemeris data of LEO satellite and Beidou navigation satellite are input into the user algorithm of Multiple Hypothesis Solution Separation (MHSS) of ARAIM, the satellite position information is output according to the satellite almanac, and then the user grid point position information is read, and the visible stars are compared with the shielding angle standard. According to the preset error model and fault mode, the fully visible satellite solution and the subset solution of each fault mode are calculated. After the threshold test of separation, the horizontal/vertical protection level and effective monitoring threshold are calculated, so as to evaluate the state of the constellation. In the basic MHSS ARAIM protection level calculation, the vertical protection value (VPL) less than 35 m is the key index to evaluate the usability, and the VPL value of this grid point is calculated by the maximum function. 
     According to the embodiment of the present invention, under the condition that the integrity risk and continuity risk are equally distributed, the vertical protection level after traversing all subset solutions and fault modes is expressed as: 
     VPL q =max ((VPL 0 ) q ,max ((VPL i ) q ), where VPL 0  is a vertical protection level under a fault-free condition, (VPL 0 ) q =K MD,q   0 ·σ 0,q +Σ n=1   Nsat |S q,n   0 |·b int,n , K MD,q   0  is an integrity and continuity risk value of the fault-free mode under a fully visible satellite subset, S q,n   0 , is a projection matrix, and b int,n  is a maximum nominal deviation of a n th  satellite; 
     VPL, is a vertical protection level corresponding to a measurement deviation of an i th  fault mode where a maximum deviation is not exceeded, (VPL i ) q =K MD,q   i ·σ i,q +Σ n=1   Nsat |S q,n   i |·b int,n +D i,q , D i,q  is a detection threshold; 
     wherein , σ 0,q =√{square root over ((GW INT G T ) q,q   −1 )}, σ i,q =√{square root over ((G T M i W INT G T ) q,q   −1 )}, G is a geometric matrix in a pseudo-range observation equation, W INT  is a parameter of a fixed error hypothesis model related to integrity, M i  is an identity matrix of a size N sat *N sat , N sat  is the number of visible satellites, and q is a q th  sample point; 
     the geometric matrix G in the pseudo-range observation equation contains the LEO satellite constellation configuration parameters x 1 , x 2 , x 3 , x 4 , in which the parameter x 1  is the orbit inclination, the parameter x 2  is the orbit altitude, the parameter x 3  is the starting value parameter of ascending intersection right ascension, and the parameter x 4  is the starting value of mean anomaly. 
     According to the present invention, the constellation configuration of the LEO augmentation system is optimized, and finally the maximum VPL is optimized to achieve the purpose of reducing the protection level. Because the protection level is optimized to the maximum value, the protection level is optimized by reducing the integrity and continuity risk value or optimizing the observation matrix and error parameters. 
     In the optimization process, the vertical protection level VPL is a function of the LEO constellation configuration parameter x 1 , x 2 , x 3 , x 4 , that is, the objective function in step 2 below. Therefore, the configuration of Low Earth Orbit LEO should be defined first. Combining with the actual launch cost, this invention selects the Walker configuration in two dimensional Lattice Flower Constellation (2D-LFC), which is a special case of 2D-LFC. 2D-LFC can define the orbit of a satellite with nine parameters, six of which are Kepler elements. 
     According to the present invention, the constraint conditions of LEO satellite constellation configuration parameters are determined, and a 2D-LFC constellation configuration will satisfy the following constraints: 
     
       
         
           
             
               
                 
                   [ 
                   
                     
                       
                         
                           N 
                           0 
                         
                       
                       
                         0 
                       
                     
                     
                       
                         
                           N 
                           C 
                         
                       
                       
                         
                           N 
                           SO 
                         
                       
                     
                   
                   ] 
                 
                 [ 
                 
                   
                     
                       
                         Ω 
                         ij 
                       
                     
                     
                       
                         Ω 
                         11 
                       
                     
                   
                   
                     
                       
                         M 
                         ij 
                       
                     
                     
                       
                         M 
                         11 
                       
                     
                   
                 
                 ] 
               
               = 
               
                 2 
                 ⁢ 
                 
                   π 
                   [ 
                   
                     
                       
                         i 
                       
                       
                         
                           - 
                           1 
                         
                       
                     
                     
                       
                         j 
                       
                       
                         
                           - 
                           1 
                         
                       
                     
                   
                   ] 
                 
               
             
             , 
           
         
       
     
     where N 0  is the number of orbital planes of LEO constellation, N SO  is the number of satellites on each orbital plane, N C , is a phase parameter, N C ∈[1, N 0 ]. Ωis a right ascension of an ascending intersection point, M is an average anomaly, where i represents an i th  orbital plane and j represents a j th  satellite. 
     Considering the external constraints, the four parameters should be constrained within a certain range, as shown in Table 1. 
     
       
         
           
               
             
               
                 TABLE 1 
               
             
            
               
                   
               
               
                 LEO constellation optimization parameters 
               
            
           
           
               
               
               
            
               
                   
                 Parameter 
                 Range 
               
               
                   
                   
               
               
                   
                 x 1  orbit inclination (rad) 
                 0-π/2 
               
               
                   
                 x 2  orbit height (km) 
                 800-2000  
               
               
                   
                 x 3  starting value of ascending 
                  0-0.782 
               
               
                   
                 intersection right ascension (rad) 
               
               
                   
                 x 4  starting value of mean anomaly 
                 0-0.3 
               
               
                   
                 (rad) 
               
               
                   
                   
               
            
           
         
       
     
     Step 2, determining objective functions of LEO satellite constellation configuration parameters x 1 , x 2 , x 3 , x 4 , eliminating calculated values of abnormal vertical protection levels, and screening initial populations of the parameters x 1 , x 2 , x 3 , x 4 ; 
     where the parameter x 1  is the orbit inclination, the parameter x 2  is the orbit altitude, the parameter x 3  is the starting value parameter of ascending intersection right ascension, and the parameter x 4  is the starting value of mean anomaly. 
     Determine the Objective Function 
     As this invention is a multiple objective optimization problem (MOOP), which is discussed based on pareto optimal solution, the screening of the initial population is particularly critical when using the Elitist Nondominated Sorting Genetic Algorithm (NSGA-II) method. 
     Firstly, the MOOP is transformed into a simple objective optimization problem (SOOP), and the objective functions of the four optimization parameters are listed as follows: 
     
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             1 
                           
                           ) 
                         
                       
                       ] 
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         0 
                       
                       ≤ 
                       
                         x 
                         1 
                       
                       ≤ 
                       
                         π 
                         2 
                       
                     
                   
                 
               
               , 
             
           
         
       
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             2 
                           
                           ) 
                         
                       
                       ] 
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         800 
                       
                       ≤ 
                       
                         x 
                         2 
                       
                       ≤ 
                       2000 
                     
                   
                 
               
               , 
             
           
         
       
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             3 
                           
                           ) 
                         
                       
                       ] 
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         0 
                       
                       ≤ 
                       
                         x 
                         3 
                       
                       ≤ 
                       0.782 
                     
                   
                 
               
               , 
             
           
         
       
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             4 
                           
                           ) 
                         
                       
                       ] 
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         0 
                       
                       ≤ 
                       
                         x 
                         4 
                       
                       ≤ 
                       0.3 
                     
                   
                 
               
               , 
             
           
         
       
     
     In mathematical processing, the minimization of the maximum value means min(max((VPL 0 ) q , max((VPL i ) q ))); because it is inconvenient for mathematical calculation, the formula takes two negative signs, and the original optimization objective is transformed into: 
       −min└−max((VPL 0 ) q , max((VPL i ) q ))┘.
 
     Then the objective functions corresponding to the four LEO satellite constellation configuration parameters x 1 , x 2 , x 3 , x 4  are described as: 
     
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             1 
                           
                           ) 
                         
                       
                       ] 
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         0 
                       
                       ≤ 
                       
                         x 
                         1 
                       
                       ≤ 
                       
                         π 
                         2 
                       
                     
                   
                 
               
               , 
             
           
         
       
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             3 
                           
                           ) 
                         
                       
                       ] 
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         800 
                       
                       ≤ 
                       
                         x 
                         2 
                       
                       ≤ 
                       2000 
                     
                   
                 
               
               , 
             
           
         
       
       
         
           
             { 
             
               
                 
                   
                     
                       min 
                       [ 
                       
                         - 
                         
                           
                             F 
                             VPL 
                           
                           ( 
                           
                             x 
                             3 
                           
                           ) 
                         
                       
                       ] 
                     
                   
                 
                 
                   
                     
                       
                         s 
                         . 
                         t 
                         . 
                             
                         0 
                       
                       ≤ 
                       
                         x 
                         3 
                       
                       ≤ 
                       0.782 
                     
                   
                 
               
               , 
             
           
         
       
       
         
           
             { 
             
               
                 
                   
                     min 
                     [ 
                     
                       - 
                       
                         
                           F 
                           VPL 
                         
                         ( 
                         
                           x 
                           4 
                         
                         ) 
                       
                     
                     ] 
                   
                 
               
               
                 
                   
                     
                       s 
                       . 
                       t 
                       . 
                           
                       0 
                     
                     ≤ 
                     
                       x 
                       4 
                     
                     ≤ 
                     0.3 
                   
                 
               
             
           
         
       
     
     where F VPL (x 1 ) is a function where a vertical protection level VPL is the parameter x 1 , F VPL (x 2 ) is a function where a vertical protection level VPL is the parameter x 2 , F VPL (X  3 ) is a function where a vertical protection level VPL is the parameter x 3 , and F VPL (x 4 ) is a function where a vertical protection level VPL is the parameter x 4 . 
     Eliminate the Calculated Values of Abnormal Vertical Protection Level 
     Due to satellite interruption, ephemeris calculation error, false alarm and other reasons, it is possible that the calculated value of individual VPL is quite different from the real value, and such VPL value cannot normally participate in the optimization calculation. Therefore, the first task is to eliminate these outliers. 
     According to the embodiment of the present invention, as shown in the initial abnormal verification diagram of the vertical protection level in  FIG.  3   , the calculated value of the abnormal vertical protection level is eliminated by the following method: 
     According to VPL, initial anomaly detection thresholds T first   + , T first   −  are set, and T first   + =μ all +2σ all , T first   − =μ all −2σ all  are defined, wherein, 
     μ all  is an average value of VPL calculated values of vertical protection levels corresponding to all sample points, and σ all  is a standard deviation of vertical protection and VPL calculated values corresponding to all sample points, 
     the calculated value of the vertical protection level VPL of each sample point is compared with T first   + , T first   − , and if the inequality T first   + ≤VPL≤T first   −  is satisfied, the sample point passes the threshold detection and waits for the initial population screening; if the inequality T first   + ≤VPL≤T first   −  is not satisfied, the sample point fails the threshold detection and is stored in a abnormal data module. 
     The failed data will be verified separately after the whole data flow, traversing all abnormal reasons, such as constellation width fault, clock ephemeris fault, signal distortion, antenna offset error, etc. If there is no match, it will be added to the initial population. 
     Screening Initial Population 
     Because the selection of initial population is the basis of genetic algorithm, the selection strategy of initial population is very important. The original random construction mode is based on a huge sample size, but the configuration parameters x 1 , x 3 , x 4  of the LEO constellation of the present invention only need to be accurate to two or three decimal places, and x 2  only needs to be accurate to one digit, otherwise, many unnecessary calculation loads will be increased. In view of this background, the present invention proposes a data initial population screening strategy based on sampling aiming at the above four parameters. 
     According to the embodiment of the present invention, the initial population screening is carried out by the following method: 
     for an orbit inclination parameter x 1 , a sampling interval is set to 0.01, and a total of 158 sample points are generated; a set of data is taken every 9 sample points as a sample of the initial population, and the 15 samples constitute the initial population of the orbit inclination parameter x 1 ; 
     for an orbit height parameter x 2 , a sampling interval is set to 1, and a total of 1200 sample points are generated; a set of data is taken every 19 sample points as a sample of the initial population, and the 60 samples form the initial population of the orbit inclination parameter x 2 ; 
     for an initial value parameter x 3  of ascending intersection right ascension, a sampling interval is set to 0.001, and a total of 79 sample points are generated; a set of data is taken every four sample points as a sample of the initial population, and the 15 samples constitute the initial population of the orbit inclination parameter x 3 ; 
     for an initial value parameter x 4  of mean anomaly, a sampling interval is set to 0.01, and a total of 30 sample points are generated; a set of data is taken every four sample points as a sample of the initial population, and the six samples constitute the initial population of the orbit inclination parameter x 4 . 
     Step 3, calculating fitness of the objective functions of the LEO satellite constellation configuration parameters x 1 , x 2 , x 3 , x 4 ; 
     After screening and sampling, an initial population with the size of N is generated, and then the rapid non-dominated sorting, selection, crossover, mutation and other operations are performed. The implementation of the algorithm is a process of survival of the fittest similar to the evolutionary theory, and the whole process is also a process of evaluating individual fitness. 
     In genetic algorithm, the fitness of individuals is used to evaluate the merits of individuals, so fitness function is needed to participate in the evaluation. Appropriate selection method of evolutionary individuals is helpful to improve the efficiency of population evolution. 
     The method calculates the fitness of an objective function of the LEO satellite constellation parameters x 1 , x 2 , x 3 , x 4 : 
     
       
         
           
             
               Fitness 
               ( 
               
                 
                   F 
                   VPL 
                 
                 ( 
                 x 
                 ) 
               
               ) 
             
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           
                             
                               F 
                               VPL 
                             
                             ( 
                             x 
                             ) 
                           
                           + 
                           
                             C 
                             min 
                           
                         
                         , 
                       
                     
                     
                       
                         
                           
                             C 
                             min 
                           
                           + 
                           
                             
                               F 
                               VPL 
                             
                             ( 
                             x 
                             ) 
                           
                         
                         &gt; 
                         0 
                       
                     
                   
                   
                     
                       
                         0 
                         , 
                       
                     
                     
                       
                         
                           
                             C 
                             min 
                           
                           + 
                           
                             
                               F 
                               VPL 
                             
                             ( 
                             x 
                             ) 
                           
                         
                         ≤ 
                         0 
                       
                     
                   
                 
                 , 
               
             
           
         
       
     
     where C min  is a preset number, and the minimum function value of the objective function F VPL (x) estimated so far is taken as the objective function, F VPL (x) is an objective function, representing a function where the vertical protection level VPL is the parameter x 1 , x 2 , x 3  or x 4 . 
     As the value of the objective function F VPL (x) meets the requirements only under the requirements of LPV-200, besides the above fitness function, the objective function should also meet the following conditions: 
     the value of the objective function F VPL (x)  F     VPL     (x)  is ≤35 m. 
     Before checking the fitness function, the requirements for the objective function should be added, and the fitness function can be calculated only after passing the limit value. 
     Step 4: after screening out an initial population of the parameters x 1 , x 2 , x 3 , x 4 , starting from a second generation population, merging a parent population and an offspring population to form a new offspring population, and performing fast non-dominated sorting, calculating crowding degrees of individuals in each non-dominated layer, randomly paring the individuals, and performing a crossover operation of a genetic algorithm between the paired two individuals. 
     According to the embodiment of the present invention, starting from the second generation, the parent population of the previous generation will be merged with the offspring population. As the sample sizes and population numbers of the four parameters are different, the present invention defines the merging ratio of the parent population and the offspring population to form a new offspring population as: 
     
       
         
           
             
               
                 pro 
                 i 
               
               = 
               
                 
                   
                     ❘ 
                     &#34;\[LeftBracketingBar]&#34; 
                   
                   
                     
                       
                         X 
                         ul 
                       
                       - 
                       
                         X 
                         ll 
                       
                     
                     
                       Δ 
                       ⁢ 
                       
                         π 
                         i 
                       
                     
                   
                   
                     ❘ 
                     &#34;\[RightBracketingBar]&#34; 
                   
                 
                 
                   2 
                   ⁢ 
                   
                     ( 
                     
                       
                         N 
                         interval 
                       
                       + 
                       1 
                     
                     ) 
                   
                 
               
             
             , 
             
               i 
               = 
               1 
             
             , 
             2 
             , 
             3 
             , 
             4 
             , 
           
         
       
     
     where X ul  is an upper limit of an optimization parameter range, X ll  is a lower limit of the optimization parameter range, Δτ i  is an sampling interval of each parameter, and N interval  is a sample interval when the initial population is generated. 
     According to different proportions, and then fast non-dominated sorting is carried out. At the same time, the crowding degree of individuals in each non-dominated layer is calculated, individuals are randomly paired, and cross operations are carried out between the paired two individuals. Usually, the crossover operator and the mutation operator cooperate with each other, and the mutation operator has strong local search ability, which makes the genetic algorithm have both global search ability and local search ability. 
     5, after the parent population and the offspring population are merged to form the new offspring population, implementing optimal preservation strategy and local optimal selection for the new offspring population, and selecting a maximum value of the objective functions as an optimal offspring, and repeating step 4 until a genetic algebra is less than a maximum genetic algebra. 
     In the evolution process of genetic algorithm, only individuals with high fitness have a chance to pass on to the next generation, while individuals with low fitness have a smaller probability of passing on to the next generation. This process of survival of the fittest is realized by selecting operators. Combined with the elite strategy, the excellent individual in the parent can enter the offspring and continue to inherit, so as to prevent the loss of Pareto optimal solution. 
     The optimal preservation strategy commonly used in the prior art can make the individuals with the highest fitness not participate in crossover and mutation, and use it to replace the individuals with the lowest fitness after crossover and mutation. By this method, the individuals with the highest fitness are retained, but this method is not easy to eliminate the local optimal solution of the algorithm, which reduces the global search ability. 
     In the initial population selection, the sampling interval Δτ i  has been used for sampling, which avoids the local optimum to a certain extent. On this basis, the local mean A m   n  of the sample is introduced for the local optimum selection, where 
     
       
         
           
             
               
                 A 
                 m 
                 n 
               
               = 
               
                 
                   
                     
                       ∑ 
                       
                         
                           5 
                           ⁢ 
                           n 
                         
                         + 
                         1 
                       
                       
                         5 
                         ⁢ 
                         
                           ( 
                           
                             n 
                             + 
                             1 
                           
                           ) 
                         
                       
                     
                     
                       x 
                       i 
                     
                   
                   k 
                 
                 ⁢ 
                 
                   ( 
                   
                     
                       n 
                       = 
                       0 
                     
                     , 
                     1 
                     , 
                     2 
                     , 
                     3 
                     , 
                     … 
                        
                     , 
                     
                       m 
                       = 
                       1 
                     
                     , 
                     2 
                     , 
                     3 
                     , 
                     … 
                        
                     , 
                     
                       [ 
                       
                         
                           [ 
                           
                             
                               
                                 X 
                                 ul 
                               
                               - 
                               
                                 X 
                                 ll 
                               
                             
                             
                               Δ 
                               ⁢ 
                               
                                 π 
                                 i 
                               
                             
                           
                           ] 
                         
                         5 
                       
                       ] 
                     
                   
                   ) 
                 
               
             
             , 
           
         
       
     
     where m represents the m th  local mean value, the maximum value is 
     
       
         
           
             
               [ 
               
                 
                   [ 
                   
                     
                       
                         X 
                         ul 
                       
                       - 
                       
                         X 
                         ll 
                       
                     
                     
                       Δ 
                       ⁢ 
                       
                         π 
                         i 
                       
                     
                   
                   ] 
                 
                 5 
               
               ] 
             
             , 
           
         
       
     
     and x i  represents the i th  sample of the corresponding parameter, it is stipulated that the average value of every five samples is taken as the local mean value, X ul  is an upper limit of the optimization parameter range, X ll  is a lower limit of the optimization parameter range, and Δτ i  is a X ul  sampling interval of each parameter. 
     According to the present invention, a selection threshold T sel  is set (determined according to the statistical average of the ARAIM protection level in practical application), after the optimal preservation strategy of the new offspring population is selected, a difference between the objective function and each local mean A m   n  is calculated, 
     if the difference is greater than the selection threshold T sel , the local optimal test is passed, and the maximum value of the objective function is taken as the optimal offspring; if the difference is less than the threshold T sel , the maximum value near the objective function will be searched as the optimal offspring. 
     After the maximum value of the objective function is selected as the optimal offspring through local optimal selection, repeat step 4 and cycle until the genetic algebra meets the end condition (less than the maximum algebra, which can be adjusted at any time according to the specific simulation environment). At this time, the vertical protection level reaches the minimum value, and the values of the corresponding optimization parameters x 1 , x 2 , x 3 , x 4  are the constellation configuration that makes the protection level of LEO satellite augmentation system the lowest. Using this configuration can ensure the availability while taking into account the economic factors. 
     According to the present situation of LEO navigation enhancement, the present invention provides a constellation configuration optimization method suitable for LEO satellite navigation enhancement system, and a constellation configuration optimization method based on ARAIM protection level algorithm is proposed by combining the non-dominated sorting genetic algorithm with elite strategy. The ARAIM protection level algorithm optimization is realized from another angle through operations such as abnormal data elimination and parameter sampling, which makes up for the integrity monitoring vacancy of LEO satellite navigation enhancement system and provides reference for future LEO satellite navigation system design and networking. 
     According to the method, the parameters to be optimized are determined by the ARAIM protection level formula, and the specific observation matrix to be optimized is obtained after the risk values are evenly distributed; Define the constraint range of the optimization parameters in the observation matrix, and optimize the objective function within the constraint range; Before optimization, abnormal data detection and data sampling are carried out to screen the initial population, and when the parent and offspring populations are merged, a specific proportion is adopted to merge them purposefully; Finally, in order to avoid local optimization, the local mean of the sample is defined, and the optimized value and the local mean are threshold-tested, so as to achieve the goal of global optimization. 
     Combined with the description and practice of the present invention disclosed here, other embodiments of the present invention will be easy to think of and understand by those skilled in the art. The true scope and gist of the present invention are defined by the claims.