Radar imaging detection satellite constellation design method

Abstract

The invention discloses a radar imaging detection satellite constellation design method, which comprises the following steps of: decomposing a constellation into a plurality of orbit surfaces, uniformly describing each orbit surface through four parameters, constructing an asymmetric constellation structure, and establishing a multi-target optimization model for solving by taking the total cost of the constellation, the maximum revisit time and the detection performance as targets. Non-dominated quantizer sorting is introduced in the sorting process, and satellite constellation design is completed. According to the method, single-satellite-based optimization of the salted fish is converted into orbital plane-based optimization, the optimization variable dimension is remarkably reduced, the method is particularly suitable for large-scale constellation design comprising hundreds of satellites, meanwhile, the Pareto frontier is obtained through non-dominated sorting, various constellation configuration schemes emphasizing different performances are provided, and the method is suitable for large-scale constellation design. Compared with an existing constellation design method, the problems that the to-be-optimized parameters are increased along with the number of satellites, and efficiency is low are solved, the scheme output by the algorithm is more reasonable through the improved sorting mode, and the method can be used for constellation design.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of radar, and particularly relates to a radar imaging satellite constellation design method. BACKGROUND

[0002] Satellite-borne radars provide opportunities for detection missions due to their wide coverage; however, a single satellite is limited by revisit time and cannot meet the needs of continuous detection missions. Therefore, the design of satellite constellations has become an effective solution.

[0003] In recent years, low earth orbit (LEO) satellite constellations have developed rapidly due to their advantages of light weight, low cost, and short development cycle. The research on using satellite-borne synthetic aperture radars (SAR) to detect low-altitude targets highlights the potential of large-scale LEO constellations in terms of detection capabilities. Through reasonable design, a low revisit time can be achieved, and even continuous coverage of a specific area can be achieved, thereby providing necessary conditions for continuous detection missions.

[0004] The design problem of constellation is usually modeled as an optimization problem, based on the parameter information of satellites, the coverage range can be calculated, and the objective functions including revisit time, coverage time, etc. are modeled, and evolutionary algorithm is used for solving. In the literature "Multi-objective optimization design of extended Walker constellation for global coverage services. IEEE International Conference on Computer and Communications, pp. 1309-1313, 2016", the orbital elements are directly taken as optimization variables, although a suitable constellation configuration is solved, but for large constellation, a large number of satellites lead to a large number of optimization variables, the search space becomes too large, and the efficiency of the algorithm cannot be improved. Referring to the existing special constellation configuration can reduce the optimization variables, the literature "Evolutionary formulations for design of heterogeneous Earth observing constellations. IEEE Aerospace Conference, pp. 1-10, 2020" adopts multiple Walker constellations to synthesize a mixed constellation, reduces the optimization variables, and verifies the effectiveness of the asymmetric structure in convergence speed and the potential of comprehensive performance. The simulation results of the above research focus on constellations of two or three tens of satellites, but the number of satellites in a large constellation often exceeds one hundred, and the existing parameter design will lead to hundreds of parameters corresponding to a large constellation, the search space is too large, and it is difficult to solve; on the other hand, the design scheme obtained by using the existing NSGA-II for solving has redundancy, and part of the solution does not meet the actual demand. SUMMARY

[0005] To solve the above technical problems, the present application provides a radar imaging detection satellite constellation design method, which solves the problem of too many parameters and slow search speed of the existing constellation design method, and realizes fast operation of constellation design.

[0006] The technical scheme adopted by the present application is as follows:

[0007] S1, taking the orbit as the core to establish the optimization parameters, and initializing the population based on the NSGA-II algorithm, and entering iteration;

[0008] First, the constellation is decomposed into N orbital planes, and 4 parameters are designed for each orbit , and an asymmetric constellation structure is constructed.

[0009] Among them, represents the number of satellites in a single orbit plane, represents the orbit inclination, represents the orbit longitude of the ascending node, represents the initial mean anomaly of the first satellite in the orbit plane, i.e. the initial phase of the orbit.

[0010] Then set the orbit semi-major axis to a fixed altitude, eccentricity to 0, argument of perigee to 0, and the satellites in the same orbit are evenly distributed according to the mean anomaly, and share the inclination and the longitude of the ascending node, i.e. the orbital elements of each satellite in the orbit are consistent with the orbit plane, then the mean anomaly calculation expression is as follows:

[0011] (1);

[0012] wherein, represents the mean anomaly of the i-th satellite; the orbital elements include: . .

[0013] Finally, based on the NSGA-II algorithm, the population is initialized and iteration is entered.

[0014] S2, based on step S1, quickly reconstruct the position information of each satellite, and convert the position information of each satellite to the station-centered coordinate system representation;

[0015] Each individual is restored to the orbital elements of the satellite, and coordinate system conversion is performed. First, the orbital elements obtained in step S1 are converted into the geocentric inertial coordinate system representation, and then converted from the geocentric inertial coordinate system to the geocentric terrestrial coordinate system, and finally converted from the geocentric terrestrial coordinate system to the station-centered coordinate system to obtain the position of the satellite in the station-centered coordinate system.

[0016] S3, a multi-objective optimization model is established to solve the total cost of the constellation, the maximum revisit time and the radar detection performance, and a non-dominated quantization sub-sorting is introduced, and screening is performed according to the sorting result until the iteration is completed, and the satellite constellation design is completed;

[0017] First, based on the coordinates of each satellite obtained in step S2, the objective function is calculated, and the first non-dominated sorting is performed according to the NSGA-II algorithm. Then a non-dominated quantization sub-sorting method is established, each objective function is classified, the second sorting is performed, and finally the sorting result is screened, the next iteration is performed, until the iteration is completed, and the satellite constellation design is completed.

[0018] Further, the step S2 is specifically as follows:

[0019] S21. Convert the six orbital roots obtained in step S1 into a geocentric inertial coordinate system representation;

[0020] First, calculate the satellite's position on the elliptical plane. At this point, consider the elliptical plane to coincide with the equatorial plane, and the coordinate axes to point in the same direction as the geocentric inertial coordinate system. The resulting orbital equation is as follows:

[0021] (2);

[0022] in, Indicates the distance between the satellite and the Earth's center. Indicates the semi-major axis of the track. This indicates the eccentricity rate.

[0023] Ignoring the Z-axis information, the satellite's position on the elliptical plane is obtained as follows:

[0024] (3);

[0025] Where T represents the transpose of the matrix.

[0026] Differentiating the position formula (3) yields the satellite velocity formula. The expression is as follows:

[0027] (4);

[0028] in, This represents the gravitational constant.

[0029] Then, by using a rotation matrix, the satellite's position in the geocentric inertial coordinate system is obtained, as shown in the following expression:

[0030] (5);

[0031] S22. Based on step S21, transform the geocentric inertial coordinate system into the geocentric-ground-fixed coordinate system;

[0032] First, let's understand Greenwich Sidereal Time (GST). The calculation expression is as follows:

[0033] (6);

[0034] in, Indicates reference time GST below; This represents the Earth's angular velocity of rotation; Indicates the current moment.

[0035] Then, by using a rotation matrix, the satellite's position in the Earth-centered Earth-fixed coordinate system is obtained, as shown in the following expression:

[0036] (7);

[0037] S23, converting from the geocentric coordinate system to the station-centered coordinate system based on step S22;

[0038] First, the longitude , latitude , and height of the observation station are determined and converted into the geocentric coordinate system, with the following expression:

[0039] (8);

[0040] wherein, represents the coordinates of the observation station, represents the radius of curvature of the earth, and represent the reference ellipsoid parameters of the earth.

[0041] Finally, the position of the satellite in the station-centered coordinate system is obtained through a rotation matrix, with the following expression:

[0042] (9).

[0043] Further, the step S3 is specifically as follows:

[0044] S31, based on the coordinates of each satellite obtained in step S2, calculate the objective function and perform the first non-dominated sorting according to the algorithm;

[0045] S311, obtain the objective function by taking the total constellation cost, the maximum revisit time, and the radar detection performance as the target;

[0046] 1) Total constellation cost ;

[0047] The total constellation cost considers the one-shot multiple satellite technology, which is simplified as the single satellite construction cost and the single orbit launch cost, with the following calculation expression:

[0048] (10);

[0049] wherein, respectively represent the total number of satellites and the total number of orbits, respectively represent the single satellite construction cost and the single orbit launch cost, respectively represent the total satellite construction cost and the total launch cost.

[0050] 2) Maximum revisit time ;

[0051] The coordinates of the target point and the satellite in the station-centered coordinate system are used to calculate the downward angle , the expression is as follows:

[0052] (11);

[0053] wherein, represents the vector from the target point to the satellite, represents the height vector of the satellite relative to the target.

[0054] Set the maximum and minimum values of the downward angle Form a ring-shaped illumination area when When it meets this interval, that is, the satellite covers the target area at this moment, denoted as The revisit time is the period when the target point is not covered by the satellite, and the maximum value in a cycle is the maximum revisit time , the expression is as follows:

[0055] (12);

[0056] wherein, represents the sampling time of the th successful illumination, represents the last successful illumination time within the observation time.

[0057] 3) Radar detection performance ;

[0058] To determine whether a certain moment meets the detection condition, the included angle of the two satellites needs to be calculated. When it is found that there are two satellites that meet the illumination condition at a certain moment, that is, meets the ring-shaped illumination area, the vectors of the two satellites relative to the target and are obtained, and the dual-baseline angle of the two satellites relative to the target is calculated .

[0059] The range that the dual-baseline angle should meet is further defined. When meets the range, it is determined that the moment meets the dual-baseline angle condition, and the expression is as follows:

[0060] (13);

[0061] wherein, respectively represent the minimum value and the maximum value of the dual-baseline angle allowed, represents the minimum value of the allowed signal-to-noise ratio, represents the current signal-to-noise ratio.

[0062] When a certain moment meets the conditions of dual-baseline angle and signal-to-noise ratio, it is recorded as a successful detection, and different constraints are designed for different requirements.

[0063] The number of successful detections and the corresponding signal-to-noise ratio in a simulation cycle are recorded, and the total signal-to-noise ratio is calculated. The total signal-to-noise ratio is taken as the reciprocal and amplified times the fitness, that is, the radar detection performance , and the expression is as follows:

[0064] (14);

[0065] wherein, indicates the signal-to-noise ratio corresponding to the sampling time of the th successful detection, indicates the time of the last successful detection in the observation time, is a constant.

[0066] 4) Objective function;

[0067] The total optimization objective and constraint expression are as follows:

[0068] (15);

[0069] S312, calculate the objective function, and perform the first non-dominated sorting according to the NSGA-II algorithm;

[0070] Based on step S311, the objective function of each individual, that is, the fitness in the NSGA-II algorithm, is obtained, and the first non-dominated sorting is performed.

[0071] First, the dominance relationship between all individuals is determined by pairwise comparison. Then, sorting begins: the first round of traversal finds all solutions that are not dominated by any other individual, and they are classified into the highest level of non-dominated layer, that is, the first front; the sorted individuals are temporarily removed from the candidate set; the above process is repeated in the remaining individuals to find a new set of non-dominated solutions, forming the second front; iteration is performed until each individual in the population is classified into a specific front, and is assigned a rank Rank according to the front it is in, wherein the front with rank 1 contains the optimal non-dominated solution in the population.

[0072] S32, establish a non-dominated quantitative sub-sorting method, grade each objective function, and perform the second sorting;

[0073] Based on step S31, when the first non-dominated sorting is completed, Rank X is formed, X is 1 for the first front, and the like, and the individual fitness in each Rank is secondarily sorted in a discrete data manner.

[0074] First, each fitness is normalized, then linear quantization is carried out, the fitness is changed from a specific value to a quantized value, and the non-dominated sorting of the individual fitness after quantization is obtained. Rank X.r, when X and r are both 1, it is the first front. Finally, based on the results of the second sorting, the congestion degree is calculated inside each front. Through the quantization operation, the small differences between individuals in the objective function are ignored when judging non-dominance, and some individuals that are originally judged as non-dominant are judged as dominant.

[0075] S33, screening according to the sorting result, the next iteration is carried out until the NSGA-II algorithm iteration is completed;

[0076] After obtaining the final sorting, it is judged whether the iteration is completed, if yes, the current result is output, if not, the population is reinitialized for iteration.

[0077] The beneficial effects of the present application: the method of the present application first decomposes the constellation into multiple orbital planes, each orbital plane is described by four parameters (inclination, right ascension of ascending node, initial phase of orbit, and number of inner satellites), and an asymmetric constellation structure is constructed. Then, a multi-objective optimization model is established to solve the problem of satellite constellation design by taking the total cost of the constellation, the maximum revisit time and the detection performance as the target. In the sorting process, the non-dominated quantization sub-sorting is introduced to improve the convergence efficiency and solution set quality. The method of the present application changes the single-star-based optimization of the present application into the orbit-based optimization, which significantly reduces the dimension of the optimization variables, and is especially suitable for large constellation design containing hundreds of satellites. At the same time, the Pareto front is obtained by non-dominated sorting, and multiple constellation configuration schemes focusing on different performances are provided. Compared with the existing constellation design method, the present application not only effectively solves the problem of excessive variable space in large-scale constellation optimization, improves the convergence speed and stability, but also achieves a good balance between cost, coverage and detection performance, and can be widely applied in the fields of air target detection, regional monitoring, and earth remote sensing. BRIEF DESCRIPTION OF DRAWINGS

[0078] Figure 1 It is a flow chart of the radar imaging detection satellite constellation design method of the present application.

[0079] Figure 2 It is a constellation parameter design schematic diagram used in the embodiment of the present application.

[0080] Figure 3 It is a non-dominated quantization sub-sorting operator schematic diagram used in the embodiment of the present application.

[0081] Figure 4 It is a non-dominated quantization sub-sorting overall framework diagram used in the embodiment of the present application.

[0082] Figure 5This is a schematic diagram showing the result of the objective function iterating with the algorithm in an embodiment of the present invention.

[0083] Figure 6 This is a schematic diagram illustrating the final result visualized using STK in an embodiment of the present invention. Detailed Implementation

[0084] The method of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0085] To facilitate understanding of the technical content of this invention by those skilled in the art, the following terms are first explained:

[0086] Term 1: Six-root number of orbits

[0087] The orbital six roots are a set of six parameters that uniquely determine the position and shape of a celestial body (such as an artificial satellite) in its orbit. They include: the semi-major axis of the orbit. The distance between the apogee of an elliptical satellite orbit and the Earth's center; eccentricity. This determines the shape of the elliptical satellite orbit; orbital inclination. The angle between the satellite's orbital plane and the Earth's equatorial plane; the right ascension of the ascending node. The geocentric angle between the satellite's ascending node and the X-axis; the argument of perigee. The angular distance from the perigee to the ascending node of an elliptical satellite orbit; the angle of true anomaly. , the angular distance from the satellite's current position to its perigee.

[0088] This invention is primarily verified using simulation experiments, and all steps and conclusions have been verified correctly using Matlab 2021. For example... Figure 1 The flowchart of a radar imaging detection satellite constellation design method of the present invention is shown below, and the specific steps are as follows:

[0089] S1. Establish the parameters to be optimized with the orbit as the core, and initialize the population based on the NSGA-II algorithm to enter the iteration;

[0090] This embodiment uses a circular orbit at the same altitude, commonly used in probe constellations, and is designed as follows: Figure 2 The individual parameters of the constellations shown are first decomposed into N orbital planes, with 4 parameters designed for each orbit. To construct an asymmetric constellation structure.

[0091] in, Indicates the number of satellites in a single orbital plane. Indicates the orbital inclination angle. The right ascension (RAAN) represents the ascending node of the orbit. This represents the initial true anomaly angle of the first satellite in the orbital plane, i.e., the initial phase of the orbit. Figure 2 middle respectively represent the number of satellites in the Nth orbital plane, the orbital inclination, the right ascension of the ascending node of the orbit, and the initial phase of the orbit.

[0092] The orbital height is determined as a hyperparameter before optimization, and the argument of perigee is meaningless under a circular orbit. Then, the semi-major axis of the orbit is set as a fixed height, the eccentricity is 0, the argument of perigee is 0, and the satellites in the same orbit are uniformly distributed according to the true anomaly and share the inclination and the right ascension of the ascending node, that is, the orbital elements of each satellite in the orbit are consistent with the orbital plane, and the true anomaly calculation expression is as follows:

[0093] (1);

[0094] wherein, the true anomaly of the Nth satellite is represented by . .

[0095] Finally, based on the NSGA-II algorithm, the population is initialized according to the simulation parameters in Table 1, and iteration is entered.

[0096] Table 1

[0097]

[0098] In the initialization stage of the embodiment, 100 sets of constellation parameters are generated, that is, 100 individuals in the algorithm, to prepare for subsequent coordinate conversion and calculation of the objective function.

[0099] In the optimization based on the orbit of the embodiment, the satellites in the same orbit are regularly distributed, that is, uniformly distributed according to the true anomaly, and share the inclination and the right ascension of the ascending node, so as to reduce the length of variables. The uniformly distributed satellites in the orbit can bring relatively stable coverage performance. In the design of the embodiment, only four parameters are used to describe one orbit, and the number of variables will not increase no matter how many satellites are in the same orbit. Therefore, the total variable length is 4*N. Compared with the existing satellite-based constellation design method, the parameter amount is only O(N), and in the one-orbit multi-satellite configuration, the parameter amount will be effectively maintained and will not change dramatically with the increase of the number of satellites.

[0100] S2, based on step S1, quickly reconstruct the position information of each satellite, and convert the position information of each satellite to the station-centered coordinate system;

[0101] Each individual is reduced to the orbital elements of the satellite, and the coordinate system is converted. First, the orbital elements obtained in step S1 are converted into the geocentric inertial coordinate system, and then converted from the geocentric inertial coordinate system to the geocentric terrestrial coordinate system, and finally converted from the geocentric terrestrial coordinate system to the station-centered coordinate system to obtain the position of the satellite in the station-centered coordinate system.

[0102] The commonly used coordinate system of the satellite position includes:

[0103] The geocentric inertial coordinate system: taking the geocenter O as the coordinate origin, the X-axis points to the vernal equinox, the Y-axis is at a 90-degree angle with the X-axis in the equatorial plane, and the Z-axis satisfies the right-hand rule.

[0104] The geocentric terrestrial coordinate system: taking the geocenter O as the coordinate origin, the X-axis points to the Greenwich meridian and the intersection of the equatorial plane, the Y-axis is at a 90-degree angle with the X-axis in the equatorial plane, and the Z-axis satisfies the right-hand rule.

[0105] The station-centered coordinate system: taking the observation station (or sensor) as the center, the X-axis points to the local geographical east direction, the Y-axis points to the local geographical north direction, and the Z-axis satisfies the right-hand rule. In the design of the embodiment, the target is taken as the center.

[0106] The geodetic coordinate system: a coordinate system established based on the reference ellipsoid surface, which uses three coordinate components of longitude L, latitude B and height H to represent the position of the ground point.

[0107] S3, a multi-objective optimization model is established to solve the target of the total cost of the constellation, the maximum revisit time and the radar detection performance, and a non-dominated quantitative sub-sorting is introduced. According to the sorting result, screening is carried out until the iteration is completed, and the satellite constellation design is completed.

[0108] Firstly, based on the coordinates of each satellite obtained in step S2, the objective function is calculated, and the first non-dominated sorting is carried out according to the NSGA-II algorithm. Then a non-dominated quantitative sub-sorting method is established, each objective function is classified, the second sorting is carried out, and finally the sorting result is screened for the next iteration until the iteration is completed, and the satellite constellation design is completed.

[0109] In the embodiment, the step S2 is specifically as follows:

[0110] S21, the orbital elements obtained in step S1 are converted into the geocentric inertial coordinate system.

[0111] Firstly, the position of the satellite on the elliptical plane is calculated, at this time the elliptical plane is regarded as coinciding with the equatorial plane, the coordinate axis is directed in accordance with the geocentric inertial coordinate system, and the expression of the orbit equation of the satellite is as follows:

[0112] (2);

[0113] wherein, denotes the distance between the satellite and the Earth's center, denotes the semi-major axis of the orbit, denotes the eccentricity.

[0114] Neglecting the Z-axis information, the position of the satellite on the elliptical plane is obtained, and the expression is as follows:

[0115] (3);

[0116] wherein, T denotes the transpose of the matrix.

[0117] Differential of the position formula (3) is obtained, that is, the satellite velocity formula , and the expression is as follows:

[0118] (4);

[0119] wherein, denotes the Earth's gravitational constant.

[0120] Then, the position of the satellite in the Earth-centered inertial coordinate system is obtained through the rotation matrix, and the expression is as follows:

[0121] (5);

[0122] S22, based on step S21, converting from the Earth-centered inertial coordinate system to the Earth-centered terrestrial coordinate system;

[0123] Firstly, the Greenwich sidereal time GST is understood , and the calculation expression is as follows:

[0124] (6);

[0125] wherein, denotes the GST at the reference time , and the J2000 time (12 o'clock on January 1, 2000) is selected in the embodiment; denotes the rotation angular velocity of the Earth; denotes the current time.

[0126] Then, the position of the satellite in the Earth-centered terrestrial coordinate system is obtained through the rotation matrix, and the expression is as follows:

[0127] (7);

[0128] S23, based on step S22, converting from the Earth-centered terrestrial coordinate system to the station-centered coordinate system;

[0129] Firstly, the longitude and the latitude ​height The expression is as follows:

[0130] (8);

[0131] wherein, represents the observation station coordinate, represents the earth curvature radius, and represents the reference ellipsoid parameters of the earth.

[0132] Finally, the satellite position in the station-centered coordinate system is obtained through a rotation matrix, and the expression is as follows:

[0133] (9).

[0134] In the embodiment, the step S3 is specifically as follows:

[0135] S31, based on the coordinates of each satellite obtained in step S2, calculate the objective function, and perform the first non-dominated sorting according to the algorithm;

[0136] S311, taking the total constellation cost, the maximum revisit time and the radar detection performance as the target, the objective function is obtained;

[0137] 1) Total constellation cost ;

[0138] The total constellation cost considers the one rocket multiple satellite technology, which is simplified as the single satellite construction cost and the single orbit launch cost, and the calculation expression is as follows:

[0139] (10);

[0140] wherein, respectively represent the total number of satellites and the total number of orbits (i.e. the number of orbit planes N), respectively represent the single satellite construction cost and the single orbit launch cost, respectively represent the total satellite construction cost and the total launch cost.

[0141] 2) Maximum revisit time ;

[0142] Using the coordinates of the target point and the satellite in the station-centered coordinate system, the downward angle is calculated, and the expression is as follows:

[0143] (11);

[0144] wherein, represents the vector from the target point to the satellite, The height vector of the satellite relative to the target.

[0145] Setting the maximum and minimum values of the lower view angle Forming a ring-shaped illumination area, when When the interval is met, i.e. the satellite covers the target area at that moment, denoted as The revisit time is the period when the target point is not covered by the satellite, and the maximum value in a cycle is the maximum revisit time , the expression is as follows:

[0146] (12);

[0147] Among them, denotes the sampling time of the th successful illumination, denotes the last successful illumination time within the observation time.

[0148] 3) Radar detection performance ;

[0149] To determine whether a certain moment meets the detection condition, the included angle of the two satellites needs to be calculated. When it is found that there are two satellites that meet the illumination condition at a certain moment, i.e. meet the ring-shaped illumination area, obtain the vectors of the two satellites relative to the target and , calculate the dual baseline angle of the two satellites relative to the target .

[0150] Further regulate the range that the dual baseline angle should meet. When meets the range, it is judged that the moment meets the dual baseline angle condition, and the expression is as follows:

[0151] (13);

[0152] Among them, denote the minimum and maximum values of the dual baseline angle respectively, denote the minimum value of the allowed signal-to-noise ratio, denote the current signal-to-noise ratio.

[0153] When a certain moment meets the conditions of dual baseline angle and signal-to-noise ratio, it is recorded as a successful detection, and different constraints are designed for different requirements.

[0154] Record the number of successful detections and the corresponding signal-to-noise ratio in a simulation cycle, and calculate the total signal-to-noise ratio. Take the reciprocal of the total signal-to-noise ratio and amplify it by times to be recorded as fitness, i.e. radar detection performance , the expression is as follows:

[0155] (14);

[0156] wherein, indicates the signal-to-noise ratio corresponding to the sampling moment of the first successful detection, indicates the moment of the last successful detection within the observation time, is a constant, and in this embodiment, 10000.

[0157] 4) the objective function;

[0158] The overall optimization objective and constraint expression are as follows:

[0159] (15);

[0160] S312, calculate the objective function, and perform the first non-dominated sorting according to the NSGA-II algorithm;

[0161] Based on step S311, the objective function of each individual, i.e., the fitness in the NSGA-II algorithm, is obtained, and the first non-dominated sorting is performed.

[0162] First, the dominance relationship between all individuals is determined by pairwise comparison. Then start sorting: the first round of traversal, find all solutions that are not dominated by any other individual, and classify them as the highest level of non-dominated layer, i.e., the first front; temporarily remove the sorted individuals from the candidate set; repeat the above process in the remaining individuals to find a new set of non-dominated solutions, forming the second front; iterate until each individual in the population is classified into a specific front, and is assigned a rank Rank according to the front it belongs to, where the front with rank 1 contains the best non-dominated solution in the population.

[0163] S32, establish a non-dominated quantitative sub-sorting method, classify each objective function, and perform the second sorting;

[0164] The non-dominated quantitative sub-sorting operator is shown in Figure 3 Based on step S31, when the first non-dominated sorting is completed, Rank X is formed, X is 1 when the first front is formed, and so on. Then, the fitness of each individual in each Rank is sorted again using the discrete data method.

[0165] First, normalize each fitness, then linearly quantize it, and the fitness is changed from a specific numerical value to a quantized numerical value. Then, the non-dominated sorting of the quantized individual fitness is performed to obtain Rank X.r, as shown in Figure 4 When X and r are both 1, it is the first front. Finally, based on the results of the second sorting, the crowding degree is calculated within each front. Through the quantization operation, the small differences in the objective function between individuals are ignored when judging non-dominance, and some individuals that were originally judged to be non-dominated are judged to be dominated.

[0166] S33, screening according to the sorting result, performing the next iteration until the iteration of the NSGA-II algorithm is completed;

[0167] After obtaining the final sorting, it is judged whether the iteration is completed, if yes, the current result is output, if not, the population is reinitialized for iteration.

[0168] The iteration process of the NSGA-II algorithm is shown in Figure 5 The line in the figure represents the average level of the individual ranked as Rank 1 in each generation. As can be seen from Figure 5 , under the constraint of the objective function, the Pareto front can still be stably reduced, and has good convergence within 100 generations.

[0169] Through MATLAB and STK joint debugging, constellation design visualization can be realized, and the algorithm result is further verified. The STK result is shown in Figure 6 .

[0170] In summary, the method of the application can effectively realize the optimization design of the high-performance low-orbit imaging detection constellation configuration of the key area. By using the parameterized modeling method based on the orbital plane, the optimization variables are converted from the single-star orbital elements to the orbital plane description parameters, which significantly reduces the dimension of the optimization variables in large constellation design and greatly improves the convergence efficiency of the algorithm. At the same time, the improved multi-objective evolutionary algorithm combines the non-dominated sorting and external storage mechanism to ensure the quality of the Pareto solution set under multiple objectives such as cost, revisit time and detection performance. The method of the application can provide an efficient and reliable constellation configuration scheme for regional continuous monitoring, air target detection, earth remote sensing and other tasks.

[0171] Those skilled in the art will appreciate that the embodiments described herein are presented for the purpose of helping the reader to understand the principles of the application and should be understood as not limiting the scope of protection of the application to such specific statements and embodiments. The application can have various modifications and changes for those skilled in the art. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the application should be included in the scope of protection of the claims of the application.

Claims

1. A radar imaging satellite constellation design method, the specific steps are as follows: S1, establish the to-be-optimized parameters based on the orbit, initialize the population based on the NSGA-II algorithm, and enter iteration; First, the constellation is decomposed into N orbital planes, each of which is designed with 4 parameters , to construct an asymmetric constellation structure; wherein denotes the number of satellites in a single orbit plane, denotes the orbit inclination, denotes the right ascension of the ascending node of the orbit, denotes the initial mean anomaly of the first satellite in the orbit plane, i.e. the initial phase of the orbit; Then set the semi-major axis of the orbit For a fixed altitude, eccentricity Is 0, the argument of perigee Is 0, and the satellites in the same orbit are uniformly distributed in true anomaly, and share the inclination and longitude of the ascending node, that is, the orbital elements of each satellite in the orbit are consistent with the orbital plane. The calculation expression of the true anomaly is as follows: （1）； wherein, denotes the true anomaly of the nth satellite; the orbital elements include: ; Finally, initialize the population based on the NSGA-II algorithm, and enter iteration; S2, quickly reconstruct the position information of each satellite based on step S1, and convert the position information of each satellite to the station-centered coordinate system; Convert each individual to the orbital elements of the satellite, and perform coordinate system conversion; first, convert the orbital elements obtained in step S1 to the geocentric inertial coordinate system, then convert from the geocentric inertial coordinate system to the geocentric terrestrial coordinate system, and finally convert from the geocentric terrestrial coordinate system to the station-centered coordinate system to obtain the position of the satellite in the station-centered coordinate system; S3, establish a multi-objective optimization model for solving with the total constellation cost, maximum revisit time, and radar detection performance as the target, introduce non-dominated quantitative sub-sorting, and perform screening according to the sorting result until the iteration ends, and the satellite constellation design is completed; First, calculate the objective function based on the coordinates of each satellite obtained in step S2, and perform the first non-dominated sorting according to the NSGA-II algorithm; then, establish a non-dominated quantitative sub-sorting method, classify each objective function, perform the second sorting, and finally perform screening according to the sorting result to perform the next iteration until the iteration ends, and the satellite constellation design is completed.

2. The radar imaging surveillance satellite constellation design method according to claim 1, wherein, The step S2 is specifically as follows: S21, convert the orbital elements obtained in step S1 to the geocentric inertial coordinate system; First, calculate the position of the satellite in the elliptical plane, at this time, the elliptical plane is regarded as coinciding with the equatorial plane, the coordinate axis is directed in the same direction as the geocentric inertial coordinate system, and the orbital equation expression of the satellite is as follows: （2）； wherein, denotes the distance of the satellite from the earth's center, denotes the semi-major axis of the orbit, denotes the eccentricity; Then, ignore the Z-axis information to obtain the position of the satellite in the elliptical plane, and the expression is as follows: （3）； Wherein, T represents the transpose of the matrix; Differentiating the position equation (3) gives the satellite velocity equation , which is expressed as follows: （4）； wherein represents the Earth's gravitational constant; Then, through the rotation matrix, the position of the satellite in the geocentric inertial coordinate system is obtained, and the expression is as follows: （5）； S22, convert from the geocentric inertial coordinate system to the geocentric terrestrial coordinate system based on step S21; First, understand the Greenwich Sidereal Time, GST The calculation expression is as follows: （6）； wherein represents the GST at the reference time point represents the rotational angular velocity of the earth represents the current time point​ Then, through the rotation matrix, the position of the satellite in the geocentric terrestrial coordinate system is obtained, and the expression is as follows: （7）； S23, convert from the geocentric terrestrial coordinate system to the station-centered coordinate system based on step S22; First, the longitude of the observation station is determined , the latitude , the altitude , and converted into the geocentric geodetic coordinate system, expressed as follows: （8）； wherein denotes the observation station coordinate, denotes the earth curvature radius, and denotes the reference ellipsoid parameter of the earth; Finally, through the rotation matrix, the position of the satellite in the station-centered coordinate system is obtained, and the expression is as follows: （9）。 3. The radar imaging surveillance satellite constellation design method of claim 1, wherein, The step S3 is specifically as follows: S31, calculate the objective function based on the coordinates of each satellite obtained in step S2, and perform the first non-dominated sorting according to the algorithm; S311, obtain the objective function with the total constellation cost, maximum revisit time, and radar detection performance as the target; 1 ) Constellation total cost ; The total constellation cost considers the one-launch-multiple-satellite technology, is simplified as the construction cost of a single satellite and the launch cost of a single orbit, and the calculation expression is as follows: （10）； wherein, Ntotal and Ntotal represent the total number of satellites and the total number of orbits, respectively, Ctotal and Ctotal represent the total cost of satellite construction and the total cost of launch, respectively, Ctotal and Ctotal represent the total cost of satellite construction and the total cost of launch, respectively, 2) Maximum revisit time ; The lower view angle is calculated by using the coordinates of the target point and the satellite in the station-centered coordinate system , and the expression is as follows: （11）； wherein, represents a vector from the target point to the satellite, represents a height vector of the satellite relative to the target; Setting maximum and minimum values for the lower viewing angle Forming a circular irradiation zone, when When this interval is met, i.e. the satellite covers the target area at that moment, it is denoted by The revisit time is the time period during which the target point is not covered by the satellite, the maximum value within a cycle, i.e. the maximum revisit time , which is expressed as follows: （12）； wherein, denotes the sampling time of the last successful illumination within the observation time; 3) Radar detection performance ; To determine whether the detection condition is met at a certain time, the included angle of the two satellites needs to be calculated. When it is found that there are two satellites satisfying the illumination condition at a certain time, the vector of the two satellites relative to the target is obtained which meets the annular illumination area and , the double-base angle of the two satellites relative to the target is calculated ; The range that the double base angle should satisfy is defined again, and when When the range is satisfied, it is judged that the double base angle condition is satisfied at this moment, and the expression is as follows: （13）； wherein respectively denote the minimum and maximum value of the double base angle allowed, denotes the minimum value of the signal-to-noise ratio allowed, denotes the current signal-to-noise ratio; When the conditions of the double-base angle and the signal-to-noise ratio are met at a certain moment, it is recorded as a successful detection, and different constraint ranges are designed for different requirements; Record the number of successful detection and the corresponding SNR in a simulation cycle, calculate the total SNR; take the reciprocal of the total SNR and amplify times as fitness, that is, radar detection performance , the expression is as follows: （14）； wherein, denotes the signal-to-noise ratio of the sampling moment corresponding to the denotes the moment of the last successful detection within the observation time, is a constant;​ 4) Objective function; The total optimization target and constraint expression are as follows: （15）； S312, calculate the objective function, and perform the first non-dominated sorting according to the NSGA-II algorithm; Based on the step S311, the objective function of each individual, i.e. the fitness in the NSGA-II algorithm, is obtained, and the first non-dominated sorting is performed; Firstly, the dominance relationship between all individuals is determined by pairwise comparison; then the sorting is started: in the first round, all solutions that are not dominated by any other individual are found and classified into the highest level of non-dominated layer, i.e. the first front; the sorted individuals are temporarily removed from the candidate set; the above process is repeated in the remaining individuals to find a new set of non-dominated solutions, forming the second front; iteration is performed until each individual in the population is classified into a specific front, and is assigned a rank according to the front it belongs to, wherein the front with rank 1 contains the optimal non-dominated solution in the population; S32, a non-dominated quantitative sub-sorting method is established, each objective function is graded, and the second sorting is performed; Based on step S31, when the first non-dominated sorting is completed, Rank X is formed, X is 1 when the first front is formed, and the like, and then the individual fitness in each Rank is sorted for the second time in a discrete data manner; Firstly, each fitness is normalized, then linear quantization is performed, the fitness is changed from a specific numerical value to a quantized numerical value, and the non-dominated sorting of the quantized individual fitness is performed to obtain Rank X.r, when X and r are both 1, it is the first front; finally, based on the result of the second sorting, the crowding degree is calculated within each front; through the quantization operation, the small differences in the objective function between individuals are ignored in the non-dominated judgment, and some individuals that are originally judged as non-dominated are judged as dominated; S33, screening is performed according to the sorting result, the next iteration is performed, and the iteration is ended until the NSGA-II algorithm is ended; After the final sorting is obtained, it is judged whether the iteration is completed, if yes, the current result is output, and if not, the population is reinitialized for iteration.

Images