Constellation optimization design method for providing navigation service for earth-moon space

Abstract

The invention discloses a constellation optimization design method for providing navigation service for earth-moon space, and belongs to the field of aerospace constellation optimization. The implementation method comprises the following steps: designing a nominal orbit of each orbit type under a CRTBP model, and generating a reference working orbit library; a plurality of performance index models are established, and comprehensive scoring is carried out based on a weighted average method. According to the earth-moon navigation constellation scheme optimization algorithm, the number of satellites in a constellation serves as the unique input, the maximum number of deployed satellites under each orbit type is restrained, simulation starting and ending time is given, after weight distribution is conducted on indexes, all constellation combination schemes are generated through a dynamic programming method, and the optimal constellation scheme is obtained. And index calculation and comprehensive evaluation sorting are carried out on each scheme to obtain an optimal scheme. An orbit number combination under an orbit family selected by an optimal scheme is used as an optimization variable, and a continuous amplitude optimization problem is simplified into a discrete number integer programming problem. And after the orbit amplitude combination is determined, the initial phase of the orbit is optimized to obtain an optimal constellation configuration.

Description

Technical Field

[0001] The present invention relates to a constellation optimization design method for providing navigation services for cis-lunar space, which is applicable to universal constellation optimization design for providing navigation services for users on the ground, on the moon and in cis-lunar space, and belongs to the technical field of aerospace constellation optimization. Background Art

[0002] With the increasing number of deep space exploration missions, especially the growing demand for lunar exploration and development, reliable and real-time communication and navigation services are crucial for the precise implementation of space missions, whether they are manned lunar missions, unmanned lunar missions or deep space exploration missions. Existing GNSS navigation satellite systems are difficult to provide complete navigation and communication services for Earth-Moon space missions, so it is necessary to establish an Earth-Moon navigation constellation that covers the entire Earth-Moon region. Constellation configuration optimization design is the key to the design of the Earth-Moon navigation constellation. Traditional methods are often used to design Earth-Moon constellations with limited mission scenarios and composed of specific orbit types. They are often not universally applicable and are difficult to be widely used in the design of Earth-Moon constellations with multiple service objects and multiple orbit types. Therefore, it is necessary to propose a more universal Earth-Moon constellation design method that targets multiple service objects and consists of multiple orbits.

[0003] Among the current design methods for the Earth-Moon constellation, the prior art 1 uses a numerical extension method to design a constellation that takes into account the periodic orbit of the Earth-Moon system L1 point as the working orbit, which cannot be applied to the design of an Earth-Moon constellation that includes other orbit types.

[0004] Prior art 2 uses a mixed integer genetic algorithm to solve the lunar communication constellation design problem of LPOs based on orbital stability and synchronization constraints, but it cannot be applied to the constellation design problems of other service objects, and the constellation evaluation indicators considered are relatively simple.

[0005] Prior art 3 uses a hypervolume ant colony optimization algorithm to solve a multi-objective optimization problem to solve the design problem of a lunar frozen orbit lunar communication and navigation constellation, but this method is also not applicable to more general constellation design problems, and it is difficult to ensure convergence for other types of constellation design problems.

[0006] In summary, the existing technology uses fewer orbit types when designing constellations, and does not comprehensively consider multiple orbit types and their combinations to analyze the navigation performance of the target area. It often only uses coverage as an evaluation indicator, and the evaluation system is relatively simple. At the same time, the above-mentioned constellation design method is highly targeted and is only effective for constellation design problems in specific scenarios. The feasibility of extending the method to constellation design in different scenarios is unknown and lacks universality. Summary of the invention

[0007] The purpose of the present invention is to solve the problems of small number of orbit types involved in constellation design, single constellation performance evaluation index and poor universality of constellation design method in the prior art, and provide a constellation optimization design method for providing navigation services for Earth-Moon space. The method first generates 17 main three-body orbit families in Earth-Moon space and the constellation reference working orbits of each family based on the three-body orbit design theory, and generates corresponding constellation working orbit libraries and reference working orbit libraries based on this. The constellation service target area is defined and gridded by cubic space or spherical space. The Earth-Moon shadow, coverage and navigation precision factor are used as constellation design constraints, and the calculation model of each constraint index is constructed, and the various indicators in the simulation scenario are statistically processed, and the constellation comprehensive evaluation model is established based on weighted average. The reference working orbit library is used in combination with the dynamic programming algorithm and the traversal method to optimize the constellation scheme according to the comprehensive evaluation model. After the constellation scheme is determined, the amplitude-phase double-layer optimization strategy based on the working orbit library is further adopted, and the heuristic optimization method such as genetic algorithm is used to solve the constellation with the best performance in the current scenario. The method provides an effective tool for the construction of the Earth-Moon space navigation system, and the optimization result provides a reliable reference for the deployment of the Earth-Moon space constellation.

[0008] The purpose of the present invention is achieved through the following technical solutions.

[0009] The present invention discloses a constellation optimization design method for providing navigation services for the earth-moon space, which is an earth-moon constellation design method based on working orbit library generation and multi-index modeling. It goes through two processes of constellation scheme determination and constellation optimization to give the optimal constellation for a specific target. First, a constellation working orbit library is established. The alternative working orbit includes 17 types of periodic orbits near DRO and five libration points. The nominal orbit of each orbit type is designed under the CRTBP model to generate a reference working orbit library. A variety of performance index models are established, including earth-moon shadow, coverage, and DOP value performance indicators, and then a comprehensive score is performed based on the weighted average method. The earth-moon navigation constellation scheme optimization algorithm takes the number of satellites in the constellation as the only input, constrains the maximum number of satellite deployments under each orbit type, gives the simulation start and end time, and assigns weights to each indicator. Then, all constellation combination schemes are generated using the dynamic programming method, and the indicator calculation and comprehensive evaluation ranking are performed for each scheme to obtain the optimal scheme. On the basis of the optimal scheme, the orbit number combination under the orbit family selected by the optimal scheme is first used as the optimization variable, and the continuous amplitude optimization problem is simplified to a discrete numbered integer programming problem. Here, the orbit number is equivalent to the amplitude. After determining the orbit amplitude combination, the orbit initial phase is further optimized. The present invention realizes the design of the Earth-Moon constellation for a customized target object, comprehensively considers multiple constellation performance indicators, maximizes the constellation performance, and gives the optimal constellation configuration.

[0010] The present invention discloses a constellation optimization design method for providing navigation services for Earth-Moon space, comprising the following steps:

[0011] Step 1: Generate constellation working orbit library B and reference working orbit library B r .

[0012] The dynamic model of the spacecraft under CRTBP is:

[0013]

[0014] in

[0015] Based on the circular restricted three-body theory, N is generated by three-body orbit design and numerical extension method. t The set of all orbits under each orbit type is the orbit family of that type, and the set of orbit families formed by all orbit types is the working orbit library B of the constellation.

[0016] The orbital families under each type include C1, C2, ..., C F (F=1,2,...,N t ) orbits, F represents the index of the orbit type. The orbital families under each orbital type are sorted and numbered according to the amplitude of the X-axis direction in the Earth-Moon rotation system. From the orbital family of the Fth orbital type, take the orbital family numbered C F / 2 is called the reference working orbit of the current orbit type. The collection of reference working orbits of all orbit types is called the reference working orbit library B r .

[0017] Step 2: Determine the target area for the Earth-Moon space constellation to provide navigation services, and perform grid processing to obtain a set of space points N tar .

[0018] When the target area is a spatial target area, the spatial target area is defined using a cubic space method, and the center coordinates P of the cubic space are given. c0 =(x c0 ,y c0 ,z c0 ), the coordinate set P of the target point in cubic space:

[0019] P c ={(x c ,y c ,z c )|x c =X0+n x ·dX+x c0 ,y c =Y0+n y ·dY+yc0 ,z c =Z0+n z ·dZ+z c0} (2)

[0020] in

[0021] Indicates rounding down. The cubic space is defined by giving the coordinate ranges parallel to the three axes in the Earth-Moon rotation system. The X direction defines the starting point X0 and the end point X f , the Y direction defines the starting point Y0 and the end point Y f , the Z direction defines the starting point Z0 and the end point Z f The interval in the X direction is dX, the interval in the Y direction is dY, and the interval in the Z direction is dZ.

[0022] According to formula (2), the cubic space is divided into equal parts to achieve gridding, and then we have

[0023] N tar =P c (3)

[0024] When the target area is a ground target area or a lunar target area, a spherical space method is used to define the ground target area and the lunar target area, and a coordinate set P of the spherical space target point is generated. s

[0025] P s ={(α,λ,R)α=α0+n α ·dα,λ=λ0+n λ ·dλ} (4)

[0026] in

[0027] Spherical space is defined in the Earth-Moon rotation system, and is described by longitude, latitude, and spherical radius. Longitude is based on the X-axis of the Earth-Moon rotation system, and is positive in the counterclockwise direction. The longitude range is α∈[0,360]°. Latitude is based on the ecliptic plane, and is positive at the North Pole. The latitude range is λ∈[-90,90]°. When defining the target area of ​​spherical space, the spherical radius R is given, the longitude starting point is α0, and the longitude end point is α f , the starting latitude is λ0, and the end latitude is λ f The longitude interval is dα and the latitude interval is dλ.

[0028] According to formula (4), the spherical space is meshed by equidistant division.

[0029] For the convenience of calculation, the spherical coordinates are converted to the rectangular coordinate system. At the same time, the coordinates of the center of the sphere are given as P s0 =(x s0 ,ys0 ,z s0 ), the coordinate set of the target point in the spherical space of any region converted to the rectangular coordinate system is updated as follows:

[0030] P s ′={(x s ,y s ,z s )|x s =Rcosλcosα+x s0 ,y s =Rcosλsinα+y s0 ,z s =Rsinλ+z s0} (5)

[0031] Then there is

[0032] N tar =P s ′ (6)

[0033] Step 3: Establish criteria for determining the Earth-Moon shadow and coverage, establish a navigation precision factor calculation model, and obtain the Earth-Moon shadow, coverage, and navigation precision factor based on the navigation precision factor calculation model.

[0034] Step 3.1, establish the criteria for determining the Earth-Moon shadow.

[0035] Let r b 、r s , r are the position vectors of the blocking celestial body, the sun, and the satellite respectively. s =r s -r,s=r b -r are the relative position vectors of the sun and satellite to the celestial body. s and R are the radii of the sun and the obstructing celestial body, respectively. The shadow axis refers to the straight line connecting the mass center of the sun and the mass center of the obstructing celestial body. The shadow determination angle is the angle between the line connecting the shadow cone vertex V1 and the satellite position and the shadow axis. The opening direction of the angle is consistent with the penumbra cone angle f1, recorded as θ. s1 is the vector from the center of the earth to the shadow cone vertex V1. The above angle and vector are determined by the following formula:

[0036]

[0037] The penumbra cone angle f1 refers to the half cone angle of the penumbra cone with V1 as the vertex. f1 is expressed as follows:

[0038]

[0039] The shadow determination conditions are as follows:

[0040]

[0041] When equation (9) is satisfied, the shadow state is e=1, otherwise, e=0.

[0042] Step 3.2: Establish coverage determination criteria.

[0043] The geometric coverage satisfies that the line connecting the satellite and the target object is not blocked by the celestial body. The visible area is divided into zone I and zone II according to the orientation of the celestial body and the target object relative to the satellite, and the angle δ between the celestial body, satellite and target is calculated:

[0044]

[0045] Among them, r sat is the position coordinate of the constellation satellite, r tar is the position coordinate of the target object, r c is the position vector of the obstructing object.

[0046] For objects r that satisfy the conditions of formula (10), tar Further judgment is made. According to the cone-shaped area formed by the satellite line of sight and the celestial body, area II is divided into three sub-areas II-1, II-2 and II-3. According to formula (11), r tar Area:

[0047]

[0048] Where d is the distance between the celestial body and the line connecting the observation position and the target position. The radius of the celestial body is R c .

[0049] When tar When located in zones II-2 and II-3, further determination of r tar Visible:

[0050]

[0051] Wherein, the satellite-target-celestial body angle δ′. Based on the above judgment criteria, when the satellite and the target are visible, the coverage state c=1, otherwise, c=0.

[0052] Step 3.3: Establish a navigation precision factor calculation model.

[0053] When calculating the navigation precision factor DOP, when the number of satellites participating in the positioning solution at a certain moment in the constellation is n t , the position vector r of the target object d =(x d ,y d ,z d ), the position vector of the i-th satellite is Calculate the pseudorange measurement residual matrix A:

[0054]

[0055] in i=1,2,...,n t The covariance matrix Q is calculated by the least squares method:

[0056]

[0057] σ represents the covariance.

[0058] The calculation model of navigation precision factor is as follows:

[0059]

[0060] The Earth-Moon shadow, coverage and navigation precision factor are obtained based on the calculation model of navigation precision factor.

[0061] Step 4: Given a time series T, perform statistical processing on the results of the Earth-Moon shadow, coverage, and navigation precision factor obtained in step 3, and use the weighted average method to establish a calculation model for the constellation comprehensive evaluation index.

[0062] Step 4.1, the time series T is expressed as:

[0063] T=[t1 t2 ... t M ] (16)

[0064] M represents the number of moments in the time series.

[0065] According to the shadow situation determined in step 3.1, the n×M-dimensional shadow matrix E is expressed as:

[0066]

[0067] where e ij represents the shadow state of the i-th satellite at the j-th moment. n represents the number of satellites in the constellation. The shadow matrix is ​​summed row by row to obtain the shadow vector E of the constellation. c :

[0068]

[0069] Then with E c The maximum value in is the shadow index S e ,Right now

[0070] S e =max[E c ] (19)

[0071] Step 4.2: According to the coverage determined in step 3.2, the coverage matrix C is expressed as:

[0072]

[0073] N is the set N in step 2 tar The number of elements of .

[0074] N c To satisfy the number of spatial points of the coverage multiplicity, Cover represents the coverage multiplicity threshold, then:

[0075]

[0076] Then the coverage ratio vector of the constellation is obtained as

[0077]

[0078] C r The average value is the coverage index S c ,Right now

[0079]

[0080] Step 4.3: According to the DOP value result calculated in step 3.3, the DOP value matrix D is expressed as:

[0081]

[0082] N d To satisfy the number of spatial points that cover the multiplicity, DOPT represents the DOP value threshold, then:

[0083]

[0084] The average of the spatial proportion of DOP values ​​in the target area that meet DOPT over time is the DOP value index S. d

[0085]

[0086] Step 4.4: Based on S in step 4.1 e , S in step 4.2 c and S in step 4.3 d The weighted average method is used to establish a calculation model for the comprehensive evaluation index of constellations.

[0087] Given a weight vector W = [w1 w2 w3], the internal elements represent shadow weight, coverage weight and DOP weight respectively. Define a performance vector S = [S e S c S d ], the comprehensive index of the constellation is

[0088]

[0089] Step 5: Use dynamic programming algorithm to generate all constellation combinations that meet basic constraints I c , based on the reference working track library B in step 1 r , the spatial point set N in step 2 tar , use the comprehensive indicators in step 4 to evaluate all constellation combinations, and take the constellation solution with the highest evaluation score as the optimal constellation solution Cons0.

[0090] Generate all constellation combinations that meet basic constraints using dynamic programming algorithm I c , expressed as

[0091]

[0092] a lF It represents the number of satellites deployed in the Fth orbit type under the lth constellation combination.

[0093] to I c Traverse and evaluate row by row, take the lth row and record it as A l =[a l1 a l2 ...a lNt ], indicating the lth constellation scheme, obtain A l The non-zero elements in the vector in accordance with From B r Get the orbit of the lth constellation and perform orbit integration according to the time series T. tar As the target, use formula (27) in step 4 to evaluate the lth constellation and get S l .

[0094] Completion I c After traversing by row, obtain N tar The optimal constellation solution for providing navigation services Cons0

[0095]

[0096] b means I c The row number where the maximum constellation evaluation score is located.

[0097] Step 6: Based on Cons0 obtained in step 5, construct the amplitude optimization problem and the phase optimization problem. Based on the constellation working orbit library B generated in step 1, combined with the spatial point set N given in step 2 tar , taking the comprehensive index of step 4 as the optimization target, the genetic algorithm is used to solve the amplitude optimization problem and the phase optimization problem in turn, and the amplitude optimal constellation Cons1 and the phase optimal constellation Cons2 are obtained in turn, that is, the navigation constellation optimization of the custom target area in the Earth-Moon space is realized.

[0098] Construct and solve the amplitude optimization problem PRO1:

[0099] Get the index vector of non-zero elements in Cons0 M0 represents the number of non-zero elements in Cons0, v o Indicates the index position of the oth non-zero element in Cons0. Get the orbital library of the vth orbital type from the constellation working orbital library B. The number of orbits in the vth type orbital library is C. v , the tracks are numbered from 1 to C v , let the index be f v , f v ∈[1,C v ] and f v ∈N + .

[0100] So the optimization variable is

[0101]

[0102] The upper and lower bounds are

[0103]

[0104] The optimization goal is

[0105]

[0106] Construct and solve the phase optimization problem PRO2:

[0107] The phase optimization process uses the initial phase angle combination of each orbit in the constellation as the optimization variable, so there is

[0108]

[0109] The upper and lower bounds are

[0110]

[0111] The optimization objective is the same as formula (32).

[0112] Genetic algorithm is used to solve PRO1 to obtain the constellation with the best amplitude, Cons1. Based on Cons1, genetic algorithm is used to solve PRO2 to obtain the constellation with the best phase, Cons2. This completes the optimization design of the navigation constellation in the custom target area of ​​the Earth-Moon space.

[0113] The method also includes step seven, obtaining the optimal navigation constellation Cons2 for the target area in the Earth-Moon space under the constellation scale constraint obtained in step six, which can provide high-precision and high-efficiency positioning and navigation services to spacecraft users in this area, thereby improving the mission execution efficiency of spacecraft users in the Earth-Moon space.

[0114] Beneficial effects:

[0115] 1. The present invention discloses a constellation optimization design method for providing navigation services for Earth-Moon space. By establishing a target area definition model, the service area of ​​the Earth-Moon navigation constellation can be customized, and on this basis, the navigation constellation is optimized, so that the method has a wide range of applications.

[0116] 2. The present invention discloses a constellation optimization design method for providing navigation services in the Earth-Moon space. Based on the weighted average method, a calculation model for the comprehensive evaluation index of the constellation is established by considering the three performance indicators of the Earth-Moon shadow, coverage and precision factor. The satellite resource configuration is evaluated by using this model to improve the reliability of the constellation system.

[0117] 3. The present invention discloses a constellation optimization design method for providing navigation services in the Earth-Moon space, and constructs a constellation optimal solution generation algorithm based on dynamic programming. The algorithm takes the constellation size as the main constraint and the maximum number of satellites deployed under each orbit type as the secondary constraint. It can handle not only single-orbit single-satellite deployment solutions, but also single-orbit multi-satellite deployment solutions, and has broad application potential.

[0118] 4. The present invention discloses a constellation optimization design method for providing navigation services in the Earth-Moon space, which combines a comprehensive evaluation model with a global optimization algorithm, uses a working orbit library to perform integer programming to optimize the orbit amplitude, simplifies the continuous amplitude optimization problem into a discrete numbered integer programming problem, avoids orbit design and correction during the optimization process, significantly improves the optimization efficiency, and realizes the optimal design of the navigation constellation in the custom target area of ​​the Earth-Moon space. BRIEF DESCRIPTION OF THE DRAWINGS

[0119] Figure 1 It is an overall flow chart of an implementation case of a constellation optimization design method for providing navigation services in Earth-Moon space disclosed by the present invention.

[0120] Figure 2 Schematic diagrams of 17 reference working tracks for implementation cases of the present invention.

[0121] Figure 3 It is the target point of the cubic space of the implementation case of the present invention.

[0122] Figure 4 It is a schematic diagram of the orbit of the optimal constellation solution for the implementation case of the present invention.

[0123] Figure 5 It is a space proportion variation curve of the optimal constellation whose DOP value meets the threshold value in the implementation case of the present invention.

[0124] Figure 6It is a space proportion variation curve of the optimal constellation whose coverage multiplicity meets the threshold value in the implementation case of the present invention. DETAILED DESCRIPTION

[0125] In order to better illustrate the purpose and advantages of the present invention, the invention is further described below with reference to the accompanying drawings and examples.

[0126] Embodiment 1:

[0127] With the center of the earth as the center, the target cubic space is defined as an area of ​​500,000 kilometers in the rotating system as the coverage target. The scale of the Earth-Moon constellation is selected as 4, and the maximum number of satellites deployed in each orbit type is 1. The simulation time is from 00:00:00 on January 1, 2030 to 00:00:00 on February 1, 2030, with a total duration of 31 days. This example uses the method proposed by the present invention to give the constellation solution with the best comprehensive performance in the above scenario.

[0128] like Figure 1 As shown, the present embodiment discloses a constellation optimization design method for providing navigation services for Earth-Moon space, and the specific implementation steps are as follows:

[0129] Step 1: Generate constellation working orbit library B and reference working orbit library B r .

[0130] The dynamic model of the spacecraft under CRTBP is:

[0131]

[0132] in

[0133] Based on the circular restricted three-body theory, N is generated by three-body orbit design and numerical extension method. t The set of all orbits under each orbit type is the orbit family of that type. The set of orbit families formed by all orbit types in the table below is the constellation's working orbit library B.

[0134] The orbital families under each type include C1, C2, ..., C F (F=1,2,...,N t ) orbits, F represents the index of the orbit type. The orbital families under each orbital type are sorted and numbered according to the amplitude of the X-axis direction in the Earth-Moon rotation system. From the orbital family of the Fth orbital type, take the orbital family numbered C F / 2 is called the reference working orbit of the current orbit type. The collection of reference working orbits of all orbit types is called the reference working orbit library B r .

[0135] In this embodiment, 17 orbit types are considered as candidate orbits for constellation operation, namely, N t = 17. The number of orbits contained in each orbit family of the working orbit library B is shown in Table 1.

[0136] Table 1 Number of orbits in each orbit family of the working orbit library

[0137]

[0138]

[0139] From each orbital family, take the F / 2 orbital composition reference working orbital library B r , as shown in Table 2. See the reference working track diagram Figure 2 .

[0140] Table 2 Reference work track library

[0141]

[0142] Step 2: Determine the target area for the Earth-Moon space constellation to provide navigation services, and perform grid processing to obtain a set of space points N tar .

[0143] When the target area is a spatial target area, the spatial target area is defined using a cubic space method, and the center coordinates P of the cubic space are given. c0 =(x c0 ,y c0 ,z c0 ), the coordinate set P of the target point in cubic space:

[0144] P c ={(x c ,y c ,z c )x c =X0+n x ·dX+x c0 ,y c =Y0+n y ·dY+y c0 ,z c =Z0+n z ·dZ+z c0} (36)

[0145] in

[0146] Indicates rounding down. The cubic space is defined by giving the coordinate ranges parallel to the three axes in the Earth-Moon rotation system. The X direction defines the starting point X0 and the end point X f, the Y direction defines the starting point Y0 and the end point Y f , the Z direction defines the starting point Z0 and the end point Z f The interval in the X direction is dX, the interval in the Y direction is dY, and the interval in the Z direction is dZ.

[0147] According to formula (2), the cubic space is divided into equal parts to achieve gridding, and then we have

[0148] N tar =P c (37)

[0149] This embodiment takes the space area as the coverage target, and the area is defined as the center of the earth, that is, P c0 =(0,0,0), a cubic space area of ​​100×100×1 million kilometers in the Earth-Moon rotation system, that is, X0=-50, X f =50, Y0=-50, Y f =50, Z0=-50, Z f =50. Given that the interval between the three axes is 50,000 kilometers, that is, dX = 5, dY = 5, and dZ = 5. The spatial point set generated by formula (36) has a total of 9261 elements. The distribution of the target points in the cubic space is shown in Figure 3 .

[0150] Step 3: Establish criteria for determining the Earth-Moon shadow and coverage, and establish a navigation precision factor calculation model.

[0151] Step 3.1, establish the criteria for determining the Earth-Moon shadow;

[0152] Let r b 、r s , r are the position vectors of the blocking celestial body, the sun, and the satellite respectively. s =r s -r,s=r b -r are the relative position vectors of the sun and satellite to the celestial body. s and R are the radii of the sun and the obstructing celestial body, respectively. The shadow axis refers to the straight line connecting the mass center of the sun and the mass center of the obstructing celestial body. The shadow determination angle is the angle between the line connecting the shadow cone vertex V1 and the satellite position and the shadow axis. The opening direction of the angle is consistent with the penumbra cone angle f1, recorded as θ. s1 is the vector from the center of the earth to the shadow cone vertex V1. The above angle and vector are determined by the following formula:

[0153]

[0154] The penumbra cone angle f1 refers to the half cone angle of the penumbra cone with V1 as the vertex. f1 is expressed as follows:

[0155]

[0156] The shadow determination conditions are as follows:

[0157]

[0158] When equation (9) is satisfied, the shadow state is e=1, otherwise, e=0.

[0159] Step 3.2, establish coverage determination criteria;

[0160] The geometric coverage satisfies that the line connecting the satellite and the target object is not blocked by the celestial body. The visible area is divided into zone I and zone II according to the orientation of the celestial body and the target object relative to the satellite, and the angle δ between the celestial body, satellite and target is calculated:

[0161]

[0162] Among them, r sat is the position coordinate of the constellation satellite, r tar is the position coordinate of the target object, r c is the position vector of the obstructing object.

[0163] For objects r that satisfy the conditions of formula (10), tar Further judgment is made. According to the cone-shaped area formed by the satellite line of sight and the celestial body, area II is divided into three sub-areas II-1, II-2 and II-3. According to formula (11), r tar Area:

[0164]

[0165] Where d is the distance between the celestial body and the line connecting the observation position and the target position. The radius of the celestial body is R c .

[0166] When tar When located in zones II-2 and II-3, further determination of r tar Visible:

[0167]

[0168] Wherein, the satellite-target-celestial body angle δ′. Based on the above judgment criteria, when the satellite and the target are visible, the coverage state c=1, otherwise, c=0.

[0169] Step 3.3, establish a navigation precision factor calculation model;

[0170] When calculating the navigation precision factor DOP, when the number of satellites participating in the positioning solution at a certain moment in the constellation is n t , the position vector r of the target object d =(x d ,y d ,zd ), the position vector of the i-th satellite is Calculate the pseudorange measurement residual matrix A:

[0171]

[0172] in This embodiment is designed for a four-star constellation, so n t =4. The covariance matrix Q is calculated by the least squares method:

[0173]

[0174] σ represents the covariance.

[0175] The calculation model of navigation precision factor is as follows:

[0176]

[0177] Step 4: Given a time series T, perform statistical processing on the results of the Earth-Moon shadow, coverage, and navigation precision factor obtained in step 3, and use the weighted average method to establish a calculation model for the constellation comprehensive evaluation index.

[0178] The time series T is expressed as:

[0179] T=[t1 t2 ... t M ] (47)

[0180] M represents the number of moments in the time series.

[0181] In this embodiment, the simulation start and end times are first converted into normalized time, and a time series T is generated with a normalized time step of 0.03. T has a total of 238 moments, that is, M=238.

[0182] Step 4.1: According to the shadow situation determined in step 3.1, the n×M-dimensional shadow matrix E is obtained and expressed as:

[0183]

[0184] where e ij represents the shadow state of the i-th satellite at the j-th moment. n represents the number of satellites in the constellation. This embodiment is a four-satellite constellation, so n=4. The shadow matrix is ​​summed row by row to obtain the shadow vector E of the constellation c :

[0185]

[0186] Then with E c The maximum value in is the shadow index S e ,Right now

[0187] Se =max[E c ] (50)

[0188] Step 4.2: According to the coverage determined in step 3.2, the coverage matrix C is expressed as:

[0189]

[0190] N is the set N in step 2 tar In this embodiment, according to step 2, it can be known that N=9261.

[0191] N c To satisfy the number of spatial points of the coverage multiplicity, Cover represents the coverage multiplicity threshold, then:

[0192]

[0193] In this embodiment, Cover=4.

[0194] Then the coverage ratio vector of the constellation is obtained as

[0195]

[0196] C r The average value is the coverage index S c ,Right now

[0197]

[0198] Step 4.3: According to the DOP value result calculated in step 3.3, the DOP value matrix D is expressed as:

[0199]

[0200] N d To satisfy the number of spatial points that cover the multiplicity, DOPT represents the DOP value threshold, then:

[0201]

[0202] In this embodiment, DOPT=20.

[0203] The average of the spatial proportion of DOP values ​​in the target area that meet DOPT over time is the DOP value index S. d

[0204]

[0205] Step 4.4: Based on S in step 4.1 e , S in step 4.2 c and S in step 4.3d The weighted average method is used to establish a calculation model for the comprehensive evaluation index of constellations.

[0206] Given a weight vector W = [w1 w2 w3], the internal elements represent shadow weight, coverage weight and DOP weight respectively. Define a performance vector S = [S e S c S d ], the comprehensive index of the constellation is

[0207]

[0208] In this embodiment, there is no indicator preference, so let W = [1 1 1].

[0209] Step 5: Use dynamic programming algorithm to generate all constellation combinations that meet basic constraints I c , based on the reference working track library B in step 1 r , the spatial point set N in step 2 tar , use formula (27) in step 4 to evaluate all constellation combinations, and take the constellation solution with the highest evaluation score as the optimal constellation solution Cons0.

[0210] Generate all constellation combinations that meet basic constraints using dynamic programming algorithm I c , expressed as

[0211]

[0212] a lF It represents the number of satellites deployed in the Fth orbit type under the lth constellation combination.

[0213] to I c Traverse and evaluate row by row, taking the lth row and recording it as Indicates the lth constellation solution, obtain A l The non-zero elements in the vector in accordance with From B r Get the orbit of the lth constellation and perform orbit integration according to the time series T. tar As the target, use formula (27) in step 4 to evaluate the lth constellation and get S l .

[0214] Completion I c After traversing by row, obtain N tar The optimal constellation solution for providing navigation services Cons0

[0215]

[0216] b means I cThe row number where the maximum constellation evaluation score is located.

[0217] In this embodiment, the scale of the Earth-Moon constellation is selected to be 4, and the maximum number of satellites deployed under each orbit type is 1. A total of 2380 constellation combinations are generated using the dynamic programming algorithm. The constellation is scored and traversed using the constellation performance evaluation. The optimal constellation solution for providing navigation services for the specified space target area is obtained as a constellation of L3 South Halo, L3 vertical, L4 vertical and L2 vertical configurations. Its orbital parameters are shown in Table 3, and the orbital schematic diagram is shown in Figure 4 shown.

[0218] Table 3 Orbital parameters of the optimal constellation solution

[0219]

[0220] The Earth-Moon shadow score of this combination is 100, the coverage score is 99.99, and the DOP score is 68.1497. The weight coefficients of the three indicators are all 1, and the comprehensive score is 89.3815. The evaluation scores of the optimal constellation solution are shown in Table 4.

[0221] Table 4 Evaluation scores of the optimal constellation solutions

[0222]

[0223] Step 6: Based on Cons0 obtained in step 5, construct the amplitude optimization problem and the phase optimization problem. Based on the constellation working orbit library B generated in step 1, combined with the spatial point set N given in step 2 tar , taking equation (27) in step 4 as the optimization target, the genetic algorithm is used to solve the amplitude optimization problem and the phase optimization problem in turn.

[0224] Construct and solve the amplitude optimization problem PRO1:

[0225] Get the index vector of non-zero elements in Cons0 M0 represents the number of non-zero elements in Cons0, v o Indicates the index position of the oth non-zero element in Cons0. Get the orbital library of the vth orbital type from the constellation working orbital library B. The number of orbits in the vth type orbital library is C. v , the tracks are numbered from 1 to C v , let the index be f v , f v ∈[1,C v ] and f v ∈N + .

[0226] So the optimization variable is

[0227]

[0228] The upper and lower bounds are

[0229]

[0230] The optimization goal is

[0231]

[0232] The genetic algorithm is used to solve PRO1 to obtain the amplitude optimal constellation Cons1.

[0233] In this example, when solving PRO1, the population size of the genetic algorithm is set to 5, the number of iterations is 50, the crossover rate is 0.9, the mutation rate is 50%, the selection pressure is 0.9, and the mutation probability is 50%. After solving, the optimal variable is [1017 8717 41045227], and the optimal performance is 93.3178. The orbital parameters corresponding to Cons1 are shown in Table 5.

[0234] Table 5 Orbital parameters of Cons1

[0235]

[0236] Construct and solve the phase optimization problem PRO2:

[0237] The phase optimization process uses the initial phase angle combination of each orbit in the constellation as the optimization variable, so there is

[0238]

[0239] The upper and lower bounds are

[0240]

[0241] The optimization objective is the same as formula (32).

[0242] Based on Cons1, the genetic algorithm is used to solve PRO2 to obtain the phase optimal constellation Cons2.

[0243] When solving PRO2 in this embodiment, the population size of the genetic algorithm is set to 5, the number of iterations is 50, the crossover rate is 0.9, the mutation rate is 50%, the selection pressure is 0.9, and the mutation probability is 50%. After solving, the optimal variable is [1017 8717 41045227], and the optimal performance is 93.3178. The orbital parameters corresponding to Cons1 are shown in Table 6. The change curve of the space proportion of the DOP value of the optimal constellation that meets the threshold is shown in Figure 5 , the change curve of the space proportion that meets the coverage threshold is shown in Figure 6 .

[0244] Table 6 Orbital parameters of Cons1

[0245]

[0246] This completes the optimization design of the navigation constellation in the custom target area in the Earth-Moon space.

[0247] Step 7: Based on the optimal navigation constellation Cons2 for the target area in the Earth-Moon space under the constellation scale constraint obtained in Step 6, it is possible to provide high-precision and high-efficiency positioning and navigation services to spacecraft users in this area, thereby improving the mission execution efficiency of spacecraft users in the Earth-Moon space.

[0248] The specific description above further illustrates the purpose, technical solutions and beneficial effects of the invention in detail. It should be understood that the above is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A constellation optimization design method for providing navigation services in Earth-Moon space, characterized in that: The steps include: Step 1: Generate constellation working orbit library B and reference working orbit library B r ; Step 2: Determine the target area for the Earth-Moon space constellation to provide navigation services, and perform grid processing to obtain a set of space points N tar ; Step 3: Establish the criteria for judging the shadow of the earth and the moon and the coverage, establish the calculation model of the navigation precision factor, and obtain the shadow of the earth and the moon, the coverage and the navigation precision factor according to the calculation model of the navigation precision factor; Step 4: Given a time series T, statistically process the results of the Earth-Moon shadow, coverage, and navigation precision factor obtained in step 3, and establish a calculation model for the constellation comprehensive evaluation index using the weighted average method; Step 5: Use dynamic programming algorithm to generate all constellation combinations that meet basic constraints I c , based on the reference working track library B in step 1 r , the spatial point set N in step 2 tar , use the comprehensive index of step 4 to evaluate all constellation combinations, and take the constellation solution with the highest evaluation score as the optimal constellation solution Cons0; Step 6: Based on Cons0 obtained in step 5, construct the amplitude optimization problem and the phase optimization problem; based on the constellation working orbit library B generated in step 1, combined with the spatial point set N given in step 2 tar , taking the comprehensive index of step 4 as the optimization target, the genetic algorithm is used to solve the amplitude optimization problem and the phase optimization problem in turn, and the amplitude optimal constellation Cons1 and the phase optimal constellation Cons2 are obtained in turn, that is, the navigation constellation optimization of the custom target area in the Earth-Moon space is realized.

2. The constellation optimization design method for providing navigation services in the cis-lunar space according to claim 1, characterized in that: The implementation method of step one is: Based on the circular restricted three-body theory, N is generated by three-body orbit design and numerical extension method. t The set of all orbits under each orbit type is the orbit family of that type, and the set of orbit families formed by all orbit types is the working orbit library B of the constellation; The orbital families under each type include C1, C2, ..., C F (F=1,2,...,N t ) orbits, F represents the index of the orbit type; sort and number the orbital families under each orbital type according to the amplitude of the X-axis direction in the Earth-Moon rotation system; take the orbital family numbered C from the orbital family of the Fth orbital type F / 2 is called the reference working orbit of the current orbit type; the set of reference working orbits of all orbit types is called the reference working orbit library B r .

3. The constellation optimization design method for providing navigation services for the cis-lunar space as claimed in claim 2, characterized in that: The implementation method of step 2 is: When the target area is a spatial target area, the spatial target area is defined using a cubic space method, and the center coordinates P of the cubic space are given. c0 =(x c0 ,y c0 ,z c0 ), the coordinate set P of the target point in cubic space: P c ={(x c ,y c ,z c )x c =X0+n x ·dX+x c0 ,y c =Y0+n y ·dY+y c0 ,z c =Z0+n z ·dZ+z c0 } (1) in Indicates rounding down; cubic space is defined by giving the coordinate ranges parallel to the three axes in the Earth-Moon rotation system. The X direction defines the starting point X0 and the end point X f , the Y direction defines the starting point Y0 and the end point Y f , the Z direction defines the starting point Z0 and the end point Z f ; The interval in the X direction is dX, the interval in the Y direction is dY, and the interval in the Z direction is dZ; According to formula (1), the cubic space is divided into equal parts to achieve gridding, and then we have N tar =P c (2) When the target area is a ground target area or a lunar target area, a spherical space method is used to define the ground target area and the lunar target area, and a coordinate set P of the spherical space target point is generated. s P s ={(α,λ,R)α=α0+n α ·dα,λ=λ0+n λ ·dλ} (3) in Spherical space is defined in the Earth-Moon rotation system, and is described by longitude, latitude, and spherical radius. The longitude starts from the X-axis of the Earth-Moon rotation system, is positive counterclockwise, and has a longitude range of α∈[0,360]°; the latitude starts from the ecliptic plane, is positive at the North Pole, and has a latitude range of λ∈[-90,90]°; when defining the target area of ​​spherical space, the spherical radius R is given, the longitude starting point is α0, and the longitude end point is α f , the starting latitude is λ0, and the end latitude is λ f ; The longitude interval is dα, and the latitude interval is dλ; According to formula (3), the spherical space is meshed by equidistant division; Convert the spherical coordinates to the rectangular coordinate system; at the same time, give the coordinates of the center of the sphere as P s0 =(x s0 ,y s0 ,z s0 ), the coordinate set of the target point in the spherical space of any region converted to the rectangular coordinate system is updated as follows: P s ′={(x s ,y s ,z s )|x s =Rcosλcosα+x s0 ,y s =Rcosλsinα+y s0 ,z s =Rsinλ+z s0 } (4) Then there is N tar =P s ′ (5)。 4. The constellation optimization design method for providing navigation services for the cis-lunar space as claimed in claim 3, characterized in that: The implementation method of step three is: Step 3.1, establish the criteria for determining the Earth-Moon shadow; Let r b 、r s , r are the position vectors of the blocking celestial body, the sun, and the satellite respectively; s s =r s -r,s=r b -r are the relative position vectors of the sun and satellite to the celestial body; R s and R are the radii of the sun and the obstructing celestial body respectively; the shadow axis refers to the straight line connecting the mass center of the sun and the mass center of the obstructing celestial body; the shadow determination angle is the angle between the line connecting the shadow cone vertex V1 and the satellite position and the shadow axis; the opening direction of the angle is consistent with the penumbra cone angle f1, recorded as θ; s1 is the vector pointing from the center of the earth to the shadow cone vertex V1; the above angle and vector are determined by the following formula: The penumbra cone angle f1 refers to the semi-cone angle of the penumbra cone with V1 as the vertex; f1 is expressed as follows: The shadow determination conditions are as follows: When equation (8) is satisfied, the shadow state e = 1, otherwise, e = 0; Step 3.2, establish coverage determination criteria; The geometric coverage satisfies that the line connecting the satellite and the target object is not blocked by the celestial body. The visible area is divided into zone I and zone II according to the orientation of the celestial body and the target object relative to the satellite, and the angle δ between the celestial body, satellite and target is calculated: Among them, r sat is the position coordinate of the constellation satellite, r tar is the position coordinate of the target object, r c is the position vector of the obstructing celestial body; For the geometrically visible object r that satisfies the condition (9), tar Further judgment is made; according to the cone area formed by the satellite line of sight and the celestial body, area II is divided into three sub-areas II-1, II-2 and II-3; according to formula (10), r tar Area: Where d is the distance between the celestial body and the line connecting the observation position and the target position; the radius of the celestial body is R c ; When tar When located in zones II-2 and II-3, further determination of r tar Visible: Wherein, the satellite-target-celestial body angle δ′; Based on the above judgment criteria, when the satellite and the target are visible, the coverage state c=1, otherwise, c=0; Step 3.3, establish a navigation precision factor calculation model; When calculating the navigation precision factor DOP, when the number of satellites participating in the positioning solution at a certain moment in the constellation is n t , the position vector r of the target object d =(x d ,y d ,z d ), the position vector of the i-th satellite is i=1,2,...,nt;Calculate the pseudorange measurement residual matrix A: in i=1,2,...,n t ; The covariance matrix Q is calculated by the least squares method: σ represents covariance; The calculation model of navigation precision factor is as follows: .

5. The constellation optimization design method for providing navigation services for the cis-lunar space as claimed in claim 4, characterized in that: The implementation method of step 4 is: Step 4.1, the time series T is expressed as: T=[t1 t2 ... t M ] (15) M represents the number of moments in the time series; According to the shadow situation determined in step 3.1, the n×M-dimensional shadow matrix E is expressed as: where e ij represents the shadow state of the i-th satellite at the j-th moment; n represents the number of satellites in the constellation; the shadow matrix is ​​summed row by row to obtain the shadow vector E of the constellation c : Then with E c The maximum value in is the shadow index S e ,Right now S e =max[E c ] (18) Step 4.2: According to the coverage determined in step 3.2, the coverage matrix C is expressed as: N is the set N in step 2 tar The number of elements of ; N c To satisfy the number of spatial points of the coverage multiplicity, Cover represents the coverage multiplicity threshold, then: Then the coverage ratio vector of the constellation is obtained as C r The average value is the coverage index S c ,Right now Step 4.3: According to the DOP value result calculated in step 3.3, the DOP value matrix D is expressed as: N d To satisfy the number of spatial points that cover the multiplicity, DOPT represents the DOP value threshold, then: The average of the spatial proportion of DOP values ​​in the target area that meet DOPT over time is the DOP value index S. d Step 4.4: Based on S in step 4.1 e , S in step 4.2 c and S in step 4.3 d The calculation model of the comprehensive evaluation index of constellations is established using the weighted average method; Given a weight vector W = [w1 w2 w3], the internal elements represent shadow weight, coverage weight and DOP weight respectively; define a performance vector S = [S e S c S d ], the comprehensive indicators of the constellation are: 。 6. The constellation optimization design method for providing navigation services in the cis-lunar space according to claim 5, characterized in that: Step 5 is implemented as follows: Generate all constellation combinations that meet basic constraints using dynamic programming algorithm I c , expressed as a lF represents the number of satellites deployed in the Fth orbit type under the lth constellation combination; to I c Traverse and evaluate row by row, taking the lth row and recording it as Indicates the lth constellation solution, obtain A l The non-zero elements in the vector in accordance with From B r Get the orbit of the lth constellation and perform orbit integration according to the time series T; tar As the target, use formula (26) in step 4 to evaluate the lth constellation and get S l ; Completion I c After traversing by row, obtain N tar The optimal constellation solution for providing navigation services Cons0 b means I c The row number where the maximum constellation evaluation score is located.

7. The constellation optimization design method for providing navigation services in the cis-lunar space according to claim 6, characterized in that: The implementation method of step six is: Construct and solve the amplitude optimization problem PRO1: Get the index vector of non-zero elements in Cons0 M0 represents the number of non-zero elements in Cons0, v o Indicates the index position of the oth non-zero element in Cons0; obtains the orbital library of the vth orbital type from the constellation working orbital library B, and the number of orbits in the vth type orbital library is C v , the tracks are numbered from 1 to C v , let the index be f v , f v ∈[1,C v ] and f v ∈N + ; So the optimization variable is The upper and lower bounds are The optimization goal is Construct and solve the phase optimization problem PRO2: The phase optimization process uses the initial phase angle combination of each orbit in the constellation as the optimization variable, so there is The upper and lower bounds are The optimization objective is the same as formula (31); The genetic algorithm is used to solve PRO1 to obtain the constellation Cons1 with the optimal amplitude. Based on Cons1, the genetic algorithm is used to solve PRO2 to obtain the constellation Cons2 with the optimal phase. This completes the optimization design of the navigation constellation in the custom target area in the Earth-Moon space.

8. A constellation optimization design method for providing navigation services for Earth-Moon space as claimed in claim 1, 2, 3, 4, 5, 6 or 7, characterized in that: The method also includes step seven, obtaining the optimal navigation constellation Cons2 for the target area in the Earth-Moon space under the constellation scale constraint obtained in step six, which can provide high-precision and high-efficiency positioning and navigation services to spacecraft users in this area, thereby improving the mission execution efficiency of spacecraft users in the Earth-Moon space.

Images