A spacecraft trajectory planning method for constellation multi-target traversal flyby

Abstract

The application discloses a spacecraft trajectory planning method for constellation multi-target traversal flyby, and belongs to the technical field of aerospace. The method comprises the following steps: introducing a flyby time constraint condition and an ideal flyby azimuth, constructing a task star flyby state optimization problem, and generating a relative motion trajectory for flying by a single target. A tangential velocity pulse is applied to an apogee of the task star to adjust a semi-major axis, and a difference in right ascension drift rate of the task star and target stars under J2 perturbation at different orbit semi-major axes is utilized to realize transfer of the task star among orbit planes of the target constellation. A correction maneuver is applied during orbit plane drift of the task star to realize matching of an amplitude angle of the task star apogee and a target star phase, so that the task star enters a resonance orbit for satellite traversal flyby in the constellation. The above-mentioned transfer trajectory optimization process of the task star among adjacent orbit planes is repeated to realize flyby traversal of a single task star to all satellites of the target constellation, and further realize spacecraft trajectory planning.

Description

TECHNICAL FIELD

[0001] The present application relates to a spacecraft trajectory planning method for constellation multi-target traversal flyby, in particular to a spacecraft trajectory planning method suitable for traversing and flying by multiple targets on the same orbit plane in a large constellation, belonging to the field of aerospace technology. BACKGROUND

[0002] Large-scale constellations have become the hardware basis for communication, navigation, remote sensing and other applications, and have gradually become a new type of space infrastructure. The spacecraft trajectory planning method for traversing multiple satellite targets in the constellation has important significance in the daily patrol and maintenance work of large-scale constellations. It not only specifies the relative motion state of the task star and the target satellite in the constellation, but also generates a patrol scheme for traversing and visiting the entire constellation, which is the key to monitoring and maintaining large-scale constellations. The current multi-target traversal technology for large-scale constellations mainly adopts the form of position and speed matching rendezvous, which has a large fuel consumption in the patrol and monitoring scenario of satellite targets in the constellation, significantly affecting the on-orbit life of the mission spacecraft. Based on this, the spacecraft trajectory planning method for constellation multi-target traversal flyby proposed in this patent uses flyby to complete the detection of targets on the same orbit plane of the constellation, and uses the Earth's non-spherical perturbation to transfer between the orbit planes of the target constellation. It not only provides a solution for the traversal and patrol of large-scale constellations, but also significantly reduces the fuel consumption of the task star compared to the rendezvous form scheme, thereby meeting the expandability requirements of more spacecraft on-orbit service tasks.

[0003] Among the developed spacecraft multi-target traversal trajectory planning methods, the first technology [1] (see: SUNG T Y, AHN J Y. Optimal deployment of satellite mega-constellation [J]. Acta Astronautica, 202 (2023): 653-669.) considers the optimization of constellation deployment, also involves the traversal strategy of multiple targets in the constellation, but uses mixed integer programming to solve the satellite deployment problem. Although it has good universality in solving general multi-target traversal problems, it has the disadvantage of large computational load in the multi-satellite flyby traversal scenario involved in this patent, making it difficult to adapt to the autonomous operation of the task star. SUMMARY

[0004] The technical problem to be solved by the spacecraft trajectory planning method for constellation multi-target traversal flyby disclosed in the application is to realize flyby of all satellites in a large-scale constellation based on Walker-d configuration in turn under the condition that fuel consumption of a mission satellite is as low as possible, to give a maneuver timing and a corresponding speed increment of the mission satellite through planning, to realize spacecraft trajectory planning for constellation multi-target traversal flyby, and to implement guidance control on the mission satellite according to the trajectory planning result of the mission satellite, so as to reduce fuel consumption required by the mission satellite for traversal and patrol between targets, and to improve execution efficiency of the patrol task without additional waiting adjustment phase.

[0005] The object of the application is achieved by the following technical scheme.

[0006] The spacecraft trajectory planning method for constellation multi-target traversal flyby disclosed in the application takes relative position and speed states of the mission satellite when the mission satellite flies by a target satellite as optimization variables, introduces flyby time constraint conditions and ideal flyby azimuth, constructs a flyby state optimization problem of the mission satellite, and generates a relative motion trajectory for flying by a single target. A tangential velocity pulse is applied at a perigee of the mission satellite to adjust a semi-major axis, and a difference in right ascension drift rate of ascending nodes of the mission satellite and the target satellite under J2 perturbation at different orbit semi-major axes is utilized to realize transfer of the mission satellite between orbit planes of the target satellite constellation. A correction maneuver is applied during orbit plane drift of the mission satellite to realize matching of the perigee amplitude angle of the mission satellite and the phase of the target satellite, so that the mission satellite enters a resonance orbit for traversal flyby of satellites in the constellation. The above transfer trajectory optimization process of the mission satellite between adjacent orbit planes is repeated to realize flyby and traversal of all satellites of the target satellite constellation by the single mission satellite, and guidance control is implemented on the mission satellite according to the trajectory planning result of the mission satellite, so as to reduce fuel consumption required by the mission satellite for traversal and patrol between targets, and to improve execution efficiency of the patrol task without additional waiting adjustment phase.

[0007] The spacecraft trajectory planning method for constellation multi-target traversal flyby disclosed in the application comprises the following steps:

[0008] Step one: The relative position and velocity state of the mission star when it flies by the target star is taken as the optimization variable. The flyby time constraint and the ideal flyby azimuth are introduced. The number of actual optimization variables is reduced by using the geometric relationship constraint of flyby time and the orbital resonance relationship between the mission star and the target star. The optimization problem of the mission star flyby state only contains 3 degrees of freedom, which converts the optimization of state vector into the optimization of space angle. The ideal flyby azimuth constraint is relaxed as the minimum objective function, and the space angle value range is established as the constraint condition to construct the mission star flyby state optimization problem. The relative position and velocity state is solved to generate the relative motion trajectory of the mission star flying by a single target.

[0009] Step 1.1 Establish the target star RTN coordinate system. According to the characteristics of the multi-target star co-orbital plane of the Walker-δ constellation configuration, the position general term expression of the target star is established by using the classical orbital elements.

[0010] The RTN coordinate system is established with the target star as the origin. The x-axis points to the outside along the radial direction of the target star, the y-axis points to the velocity direction in the target star orbital plane, and the z-axis is perpendicular to the target star orbital plane and forms a right-hand system with the x-axis and the y-axis.

[0011] In the Walker-δ constellation configuration, the orbital elements of the target stars in the same orbital plane only differ in the true anomaly. At this time, the orbital elements of the jth target satellite are represented by the vector as shown in formula (1)

[0012]

[0013] In the formula: a t represents the semi-major axis of the orbital plane, e t represents the eccentricity, i t represents the inclination, Ω t represents the eccentricity, ω t represents the argument of perigee, θ0represents the true anomaly at the initial time, and N represents the total number of satellites on this orbital plane.

[0014] Step 1.2 Take the moment when the distance between the mission star and the target star is the shortest as the flyby time t s . The relative position vector δr s and the relative velocity vector δv s of the mission star and the target star at this moment are taken as the optimization variables. According to the load capacity of the mission star and its mission purpose, the relative position and relative velocity constraints of the mission star when flying by the target star are set.

[0015] The mission star state δr s and δvs satisfy the following constraints

[0016] δr s = ||δr s || = r fix (2)

[0017] δv s = ||δv s || ∈ [v L , v U ] (3)

[0018] where: r fix is the fixed distance between the target and the mission star, v L , v U are the lower and upper bounds of the flyby relative velocity.

[0019] The ideal flyby orientation of the mission star is denoted by the azimuth angle α t and the elevation angle β t , and the corresponding unit vector is denoted by l tar . It is written as the ideal flyby orientation vector [l tar ] t as follows

[0020] [l tar ] t = [cosβ t cosα t cosβ t sina t sinβ t ] T (4)

[0021] Step 1.3. The optimization variables described in step 1.2 are reduced in dimension and converted using the constraints between the mission star and the target at the flyby time. δr s is described by the azimuth angle α and the elevation angle β as a spatial vector, and δv s is described by the direction angle γ and the relative velocity magnitude δv s as a vector perpendicular to the plane of δr s . The resonance ratio of the mission star orbit and the target orbit is determined by the number of satellites contained in each orbit plane of the target constellation, and the semi-major axis of the mission star orbit is determined. The relationship between the semi-major axis and the relative velocity δv s is constructed using the energy equation, and the number of actual optimization variables is further reduced, so that the established flyby state optimization problem of the mission star contains only 3 degrees of freedom. The ideal flyby orientation constraint established in step 1.2 is relaxed as the objective function, and the actual flyby position δr s and the ideal flyby orientation l tarThe minimum angle between them is the optimization target, and the construction of the state optimization problem of the mission star flyby is completed. The reduction of the number of optimization variables and the relaxation of the ideal flyby azimuth constraint condition can make the established mission star flyby state optimization problem easier to solve, accelerate the iteration convergence speed of the optimization solver, and improve the problem convergence.

[0022] Let t s be the time when the target star is in the inertial system, and the position and velocity state vector of the target star is denoted as r ts and v ts . Let δr s be expressed by two angles α and β in the RTN coordinate system of the target star, then the relative position [δr s ] t in the RTN coordinate system is expressed as

[0023] [δr s ] t = δr s · [cosβcosα cosβsina sinβ] T (5)

[0024] At the same time, by satisfying the minimum distance necessary condition as described in equation (6), it is ensured that the mission star and the target star are in the closest distance position

[0025] δr s ·δv s = 0 (6)

[0026] Equation (6) imposes additional constraints on δv s , and its coordinate components [δv s ] t in the RTN coordinate system are expressed as

[0027] [δv s ] t = δv s L z (-α)L y (β)[0 cosγ sinγ] T (7)

[0028] Where L x , L y , L z are rotation matrices along the x, y, z axes.

[0029] The state of the tracking star at t s should be

[0030]

[0031] In the formula: r cs , r ts respectively represent ts Positions of the task star and the target star at the moment, v cs , v ts , respectively, represent the velocities of the task star and the target star, [δr s ] eci , [δv s ] eci are the component representations of the relative position and the relative velocity vector in the earth-centered inertial system, obtained by coordinate transformation from equations (5) and (7)

[0032] [δr s ] eci = L z (-Ω t ) L x (-i t ) L z (-θ t -ω t )[δr s ] t (9)

[0033] [δv s ] eci = L z (-Ω t ) L x (-i t ) L z (-θ t -ω t )[δv s ] t (10)

[0034] To satisfy the same relative state through the j+1 target star, the orbit period T c of the task star should satisfy

[0035]

[0036] In the equation, T t represents the orbit period of the target star, and N represents the total number of satellites on the orbit plane. Further, the semi-major axis a c of the tracking star is calculated as

[0037]

[0038] In the equation, μ is the earth-centered gravitational constant. From the energy equation, equation (8) satisfies

[0039]

[0040] In the equation, r cs = ||r cs ||, v cs = ||vcs ||.

[0041] The performance index function is defined as the angle between the ideal direction vector l tar and the actual flyby direction [δr s ] t / ||δr s ||, denoted as

[0042] J = min [arccos (l tar · δr s / ||δr s ||)] (14)

[0043] Based on the above analysis and modeling process, the state optimization problem of the mission spacecraft flying by the target spacecraft has the following four free variables

[0044]

[0045] Satisfying one constraint (13), a 3-degree-of-freedom state optimization problem of the mission spacecraft flying by the target spacecraft is formed.

[0046] Step 1.4 Solve the state optimization problem of the mission spacecraft flying by the target spacecraft established in step 1.3 to obtain the relative state of the mission spacecraft when flying by the target spacecraft, and form the trajectory of the mission spacecraft flying by a single target spacecraft. Due to the introduction of the resonance orbit constraint condition, the target spacecraft can naturally realize the traversal flyby of all target spacecrafts on the same orbital plane on the orbit of flying by a single target spacecraft.

[0047] As an optimization, the fmincon function in MATLAB is used to solve the above optimization problem.

[0048] Step two: Take the time when the mission spacecraft flies by the last target in the single orbital plane of the target constellation according to the relative motion trajectory optimized in step one as the transfer departure time, adjust the semi-major axis by applying a tangential velocity pulse, and realize the orbital plane transfer of the mission spacecraft by taking advantage of the difference in the right ascension drift rate of the mission spacecraft and the target spacecraft at different orbital semi-major axes under J2 perturbation. In order to realize that the mission spacecraft is exactly at the perigee position after completing the orbital plane transfer, according to the relationship that the drift transfer time is an integer multiple of the mission spacecraft period, an equation is constructed with the tangential velocity pulse applied at the transfer departure time as the only independent variable, and the corresponding velocity pulse control quantity is solved to realize that the mission spacecraft is exactly at the perigee position after completing the orbital plane transfer. The mission spacecraft realizes the transfer between adjacent orbital planes in the target constellation by taking advantage of the difference in perturbation forces between the mission spacecraft and the satellites in the target constellation, without consuming any fuel for maneuvering, which will significantly prolong the on-orbit service life of the mission spacecraft.

[0049] The time when the mission spacecraft applies a velocity pulse at the perigee is denoted as t2, and the time when the mission spacecraft reaches the adjacent orbital plane and flies by the first target spacecraft is denoted as t3.

[0050] At time t3, the mission star should be at the perigee position of the orbit, i.e. the orbit plane drift time of the mission star is an integer multiple of the mission star period. The above integer multiple relationship is expressed as

[0051]

[0052] In the formula, the mod(·) function represents the remainder, ΔT is the orbit plane shift time obtained in formula (19), T c1 represents the orbit period of the mission star after the pulse is applied at time t1. Where ΔT is given by the difference between the ascending node precession rates of the mission star and the target star. Under the influence of J2 perturbation, the ascending node right ascension will drift according to the Gauss perturbation equation as shown in formula (17):

[0053]

[0054] In the formula, a c represents the semi-major axis of the mission star after the tangential velocity increment is applied, i c represents the inclination of the orbit of the mission star. Similarly, the orbit plane drift speed of the target star is

[0055]

[0056] From formula (17) and formula (18), by changing the semi-major axis a c of the mission star, the difference between the drift speeds of the ascending nodes of the mission star and the target star is used to complete the transfer of the mission star between adjacent target orbit planes. The difference in the right ascensions of the ascending nodes between adjacent orbit planes of the target star constellation is denoted as ΔΩ, and the time required to complete the drift is

[0057]

[0058] The tangential velocity increment applied at time t2 is denoted as Δv pe , and formula (17) to formula (19) are brought into formula (16) to obtain an equation with Δv pe as the only unknown. Solving the equation group constructed by formula (16), formula (17) to formula (19) obtains the corresponding velocity pulse control quantity, and after the orbit plane transfer is completed, the mission star is exactly at the perigee position of the orbit. Δv pe is called the orbit-raising velocity increment. The mission star uses the difference in perturbation forces between itself and the satellites in the target star constellation to achieve transfer between adjacent orbit planes in the target star constellation, without consuming any fuel for maneuvering, which will significantly prolong the on-orbit service life of the mission star.

[0059] As a preferred solution, the fsolve function in MATLAB is used to solve and obtain the maneuvering pulse that needs to be applied at time t2. ​

[0060] Step three: Based on the tangential velocity impulse control amount of the transfer departure time obtained in step two, a correction maneuver is applied in the process of the mission star orbit plane drift to achieve the matching of the mission star perigee amplitude angle and the target star phase. Step one is repeated to calculate the relative position and velocity state of the mission star when it flies by the target. An optimization problem is constructed with the minimum correction maneuver velocity increment as the objective function, the relative position of the target star that satisfies the minimum flyby perigee angular distance of the mission star terminal position as the constraint, and the time of applying the correction maneuver as the optimization variable. The relationship between the time of applying the correction maneuver and the size of the correction maneuver velocity increment is established through the Lambert problem, and the optimization problem is solved to obtain the control amount of the intermediate correction maneuver applied in the process of the mission star orbit plane transfer. Since the mission star and several satellites on the adjacent orbit plane all have matching opportunities, under the premise of satisfying the constraint that the mission star is at the perigee after the transfer is completed, only the intermediate correction needs to be applied to achieve the coincidence of the mission star and a certain target star after the orbit plane transfer is completed.

[0061] The intermediate correction velocity increment Δv inter is applied to satisfy the following terminal constraint: At t3, the mission star orbit perigee coincides with a certain target star. The terminal constraint is expressed as

[0062]

[0063] In the formula: ω c3 represents the perigee amplitude angle of the mission star at t3, u i represents the latitude amplitude angle of the i-th target star on the 2nd orbit plane, and N represents the number of target stars on the target orbit.

[0064] The target star i * with the minimum perigee angular distance to the mission star at t3 is found, i.e.

[0065]

[0066] Step one is repeated to obtain the traversal flyby mission orbit with the target star numbered i * as the starting point, which is the target mission orbit of the mission star. The perigee of the target mission orbit is the target position of the mission star at t3, denoted as

[0067] Considering that the intermediate correction time t2 is placed within the previous orbit of the mission star reaching the target mission orbit, the correction impulse Δv inter is locally optimized to minimize the impact of the intermediate correction on the orbit plane drift time, i.e.

[0068] min J' = || Δv inter || (22)

[0069] satisfy the following constraints

[0070]

[0071] t3-t2≤T c1 (24)

[0072] where: r c3 represents the target star position at t3, T c1 represents the orbit period of the target star after applying Δv pe The Gauss algorithm is used to solve the Lambert problem, and the fmincon function is combined to solve the above local optimization problem.

[0073] Step four: when the target star reaches the flyby target star relative position obtained in step three, a pulse Δv pe2 is applied to make the target star relative velocity satisfy the flyby condition, so that the target star enters the resonance orbit of the satellite flyby in the constellation, and the applied velocity pulse is recorded as the on-orbit velocity increment control quantity.

[0074] Step five: the on-orbit velocity increment obtained in step two, the correction velocity increment obtained in step three, and the on-orbit velocity increment obtained in step four are sequentially connected in time sequence to generate a control time sequence instruction for the target star to flyby all satellites in adjacent orbital planes of the target star constellation, so as to realize the spacecraft trajectory planning of the constellation multi-target flyby. According to the trajectory planning result of the target star, the guidance control of the target star can reduce the fuel consumption required for the target star to flyby and patrol between targets, and does not need to perform additional waiting adjustment phase, thereby improving the execution efficiency of the patrol task.

[0075] Advantages:

[0076] 1. The spacecraft trajectory planning method for constellation multi-target flyby disclosed by the application reduces the number of actual optimization variables by using the geometric relationship constraint at the flyby time and the orbit resonance relationship between the target star and the target star, so that the flyby state optimization problem of the target star only contains three degrees of freedom, and the flyby state optimization problem of the target star is easier to solve.

[0077] 2. The spacecraft trajectory planning method for constellation multi-target flyby disclosed by the application converts the optimization of the state vector into the optimization of the space angle, and relaxes the ideal flyby azimuth constraint as the minimum objective function, so that the flyby state optimization problem of the target star has better convergence, and the solving speed of the optimization solver can be accelerated.

[0078] 3, The spacecraft trajectory planning method for constellation multi-target traversal flyby disclosed by the application can be applied to any constellation target based on Walker-d configuration, and has no requirements for the number of orbit planes included in the constellation and the number of satellites on each orbit plane. Since the large-scale constellation currently implemented and deployed is based on Walker-d configuration, the method has certain universality.

[0079] 4, The spacecraft trajectory planning method for constellation multi-target traversal flyby disclosed by the application can make the mission star directly enter the traversal flyby mission orbit after completing the drift between adjacent orbit planes of the constellation through one-time midcourse correction, without additional waiting or phase adjustment, so that the mission execution efficiency is improved.

[0080] 5, The spacecraft trajectory planning method for constellation multi-target traversal flyby disclosed by the application realizes the transfer of the mission star between adjacent orbit planes of the target constellation through the difference in right ascension drift rate of the ascending node caused by the different orbit semi-major axes of the mission star and the target star under J2 perturbation, which can reduce the fuel consumption when adjusting the orbit plane and significantly prolong the on-orbit service life of the mission star. BRIEF DESCRIPTION OF DRAWINGS

[0081] Figure 1 is a flowchart of the spacecraft trajectory planning method for constellation multi-target traversal flyby disclosed by the application;

[0082] Figure 2 is the target star RTN coordinate system established in step one of the embodiment 1;

[0083] Figure 3 is the relative motion trajectory of the mission star for single flyby of the target star in the embodiment 1;

[0084] Figure 4 is the state of the mission star when transferred to the adjacent target mission orbit in the embodiment 1. DETAILED DESCRIPTION

[0085] In order to better illustrate the purposes and advantages of the application, the application is explained in detail below by optimizing the trajectory of the mission star in the constellation traversal flyby scenario.

[0086] As shown in Figure 1 The spacecraft trajectory planning method for constellation multi-target traversal flyby disclosed by the embodiment, the specific implementation steps are as follows:

[0087] Step one: The relative position and velocity state of the mission star when it flies by the target star is taken as the optimization variable. The flyby time constraint and the ideal flyby azimuth are introduced. The number of actual optimization variables is reduced by using the geometric relationship constraint of flyby time and the orbital resonance relationship between the mission star and the target star. The optimization problem of the flyby state of the mission star only contains 3 degrees of freedom, which converts the optimization of the state vector into the optimization of the space angle. The ideal flyby azimuth constraint is relaxed as the minimum objective function, and the space angle value range is established as the constraint condition to construct the flyby state optimization problem of the mission star, and the relative position and velocity state is obtained to generate the relative motion trajectory of the mission star flying by a single target.

[0088] Step 1.1 Establish the target star RTN coordinate system. According to the characteristics of the multi-target star co-orbital plane of the Walker-δ constellation configuration, the position general term expression of the target star is established by using the classical orbital elements.

[0089] The RTN coordinate system is established with the target star as the origin. The x-axis points to the outside along the radial direction of the target star, the y-axis points to the velocity direction in the target star orbital plane, and the z-axis is perpendicular to the target star orbital plane and forms a right-hand system with the x-axis and the y-axis.

[0090] In the Walker-δ constellation configuration, the orbital elements of the target stars in the same orbital plane only differ in the true anomaly. At this time, the orbital elements of the jth target satellite are represented by the vector as shown in equation (25)

[0091]

[0092] In the formula: a t represents the semi-major axis of the orbital plane, e t represents the eccentricity, i t represents the inclination, Ω t represents the eccentricity, ω t represents the argument of perigee, θ0 represents the true anomaly at the initial time, and N represents the total number of satellites on this orbital plane.

[0093] Step 1.2 Take the moment when the distance between the mission star and the target star is the shortest as the flyby time t s . The relative position vector δr s and the relative velocity vector δv s of the mission star and the target star at this moment are taken as the optimization variables. According to the payload capacity of the mission star and its mission purpose, the relative position and relative velocity constraints of the mission star when it flies by the target star are set.

[0094] The mission star state δr s and δvs satisfy the following constraints

[0095] δr s = ||δr s || = r fix (26)

[0096] δv s = ||δv s || ∈ [v L , v U ] (27)

[0097] where: r fix is the fixed distance between the target star and the task star, v L , v U are the lower bound and upper bound of the flyby relative velocity, respectively.

[0098] The ideal orientation of the task star when flying by the target is set, which is represented by the azimuth angle α t and the elevation angle β t , and the corresponding unit vector is denoted as l tar . It is written as the ideal flyby orientation vector of the target star RTN coordinate system [l tar ] t , as follows

[0099] [l tar ] t = [cosβ t cosα t cosβ t sina t sinβ t ] T (28)

[0100] Assuming there is a target star constellation with Walker-δ as the basic configuration, all target satellites on two adjacent orbital planes are implemented flyby traversal using the spacecraft trajectory planning method for constellation multi-target flyby of the present patent disclosure. Two orbital planes of the target star constellation are involved in this embodiment, both with an orbital height of 550 km and an orbital inclination of 53°, and the ascending node right ascension is 0° and 5°, respectively. The orbital plane with an ascending node right ascension of 0° is referred to as the first orbital plane, and the orbital plane with an ascending node right ascension of 5° is referred to as the second orbital plane. There are 22 target stars evenly distributed on each orbital plane, and the target on the second orbital plane lags behind the corresponding target on the first orbital plane by 8.8° in phase. The initial position of the task star is set to fly by the corresponding satellite on the first orbital plane, and it is set to phase 0°.

[0101] Step 1.3 uses the constraints between the task star and the target at the flyby time to reduce the dimension and convert the optimization variables described in step 1.2, and δrs As a space vector, it is described by azimuth angle a and elevation angle β, δv s Then as a vector perpendicular to δr s in the plane, it is described by the direction angle γ and the relative velocity magnitude δv s . The resonance ratio between the mission satellite orbit and the target satellite orbit is determined by the number of satellites contained in each orbit plane of the target constellation, and then the semi-major axis of the mission satellite orbit is determined. The relationship between the semi-major axis and the relative velocity δv s is constructed using the energy equation, and the number of actual optimization variables is further reduced, so that the established flyby state optimization problem of the mission satellite only contains 3 degrees of freedom. The ideal flyby azimuth constraint established in step 1.2 is relaxed as the target function, and the angle between the actual flyby position δr s and the ideal flyby azimuth l tar is minimized as the optimization goal, and the construction of the flyby state optimization problem of the mission satellite is completed. The above reduction of the number of optimization variables and relaxation of the ideal flyby azimuth constraint condition can make the established flyby state optimization problem of the mission satellite easier to solve, speed up the iteration convergence speed of the optimization solver, and improve the problem convergence.

[0102] The position and velocity state vector of the target satellite in the inertial system at time t s is denoted as r ts and v ts . The δr s is expressed by two angles a and β in the RTN coordinate system of the target satellite, and the relative position [δr s ] t in the RTN coordinate system is expressed as

[0103] [δr s ] t = δr s · [cosβcosα cosβsina sinβ] T (29)

[0104] At the same time, by satisfying the minimum distance necessary condition as shown in equation (30), it is ensured that the mission satellite and the target satellite are at the closest distance position

[0105] δr s · δv s = 0 (30)

[0106] Equation (30) imposes additional constraints on δv s , and its coordinate components [δv s ] t in the RTN coordinate system are expressed as

[0107] [δv s ] t = δv s Lz (-α)L y (β)[0 cosγ sinγ] T (31)

[0108] where L x , L y , L z are rotation matrices along x, y, z axes.

[0109] The state of the chaser at time t s should be

[0110]

[0111] where r cs , r ts are the positions of the chaser and the target at time t s , v cs , v ts are the velocities of the chaser and the target, [δr s ] eci , [δv s ] eci are the relative position and velocity vectors in the ECI frame, which are obtained from eqs. (5), (7) by coordinate transformation

[0112] [δr s ] eci = L z (-Ω t ) L x (-i t ) L z (-θ t -ω t )[δr s ] t (33)

[0113] [δv s ] eci = L z (-Ω t ) L x (-i t ) L z (-θ t -ω t )[δv s ] t (34)

[0114] In order to satisfy the same relative state through the j+1 target, the orbital period T c of the chaser should satisfy

[0115]

[0116] where T is the orbital period of the target satellite, N is the total number of satellites in the same orbital plane. Then the semi-major axis a of the tracking satellite is calculated t c

[0117]

[0118] where μ is the gravitational constant. According to the energy equation, equation (8) satisfies

[0119]

[0120] where r cs =||r cs ||, v cs =||v cs ||.

[0121] The performance index function is defined as the angle between the ideal direction vector l tar and the actual flyby direction [δr s ] t / ||δr s ||, denoted as

[0122] J = min [arccos (l tar · δr s / ||δr s ||)] (38)

[0123] Based on the above analysis, the state optimization problem of the mission satellite flying by the target satellite has the following four free variables

[0124]

[0125] satisfying one constraint equation (13), which constitutes a 3-degree-of-freedom state optimization problem of the mission satellite flying by the target satellite. Considering the concealment when flying by the target satellite, the flyby point is set in the zenith direction of the target satellite, i.e. α t = 0, β t = 0. The relative distance r fix of the mission satellite flying by the target satellite is set to 50 km, and the upper and lower bound constraints of the relative velocity v L = 100 m / s, v U = 200 m / s.

[0126] Step 1.4. Solve the state optimization problem of the mission satellite flying by the target satellite established in step 1.3 to obtain the relative state of the mission satellite when flying by the target satellite, which constitutes the trajectory of the mission satellite flying by a single target satellite. Due to the introduction of the resonance orbit constraint condition, the target satellite can naturally realize the traversal flyby of all target satellites in the same orbital plane on the orbit of flying by a single target satellite. ​​

[0127] The above optimization problem is solved using the fmincon function in MATLAB.

[0128] Step two: The mission star flies over the last target in the constellation single orbital plane at the time of the transfer departure time according to the relative motion trajectory optimized in step one, and a tangential velocity pulse is applied to adjust the semi-major axis, and the difference in the ascending node right ascension drift rate of the mission star and the target star under J2 perturbation at different orbital semi-major axes is used to realize the orbital plane transfer of the mission star. In order to realize that the mission star is exactly at the perigee position after completing the orbital plane transfer, according to the relationship that the drift transfer time of the mission star is an integer multiple of the mission star period, an equation is constructed with the tangential velocity pulse applied at the transfer departure time as the only independent variable, and the corresponding velocity pulse control quantity is obtained by solving the equation, so that the mission star is exactly at the perigee position after completing the orbital plane transfer. The mission star uses the difference in perturbation forces between itself and the satellites in the target constellation to realize the transfer between adjacent orbital planes in the target constellation, without consuming any fuel for maneuvering, which will significantly prolong the on-orbit service life of the mission star.

[0129] The time when the mission star applies a velocity pulse at the perigee is recorded as t2, and the time when the mission star reaches the adjacent orbital plane and flies over the first target star is recorded as t3.

[0130] At t3, the mission star should be at the perigee position of the orbit, that is, the orbital plane drift time of the mission star is an integer multiple of the mission star period. The above integer multiple relationship is expressed as

[0131]

[0132] In the formula: the mod(·) function represents the remainder, ΔT is the orbital plane transfer time obtained in formula (19), T c1 represents the orbital period of the mission star after applying a pulse at t1. Where ΔT is given by the difference between the ascending node progression rates of the mission star and the target star. Under the influence of J2 perturbation, the ascending node right ascension will drift according to the Gauss perturbation equation as shown in formula (17):

[0133]

[0134] In the formula: a c represents the semi-major axis of the mission star after applying a tangential velocity increment, i c represents the inclination of the mission star orbit. Similarly, the orbital plane drift speed of the target star is

[0135]

[0136] From formula (41) and formula (42), by changing the semi-major axis a c ​, the transfer of the mission satellite between adjacent orbital planes of the target constellation is completed by using the difference between the drift velocities of the mission satellite and the ascending node of the target satellite. The difference in right ascension of the ascending node between adjacent orbital planes of the target constellation is denoted as ΔΩ, and the time required to complete the drift is

[0137]

[0138] The tangential velocity increment applied at time t2 is denoted as Δv pe , and the equation with Δv pe as the only unknown is obtained by substituting equations (41) to (43) into equation (40). The equation group constructed by equations (40), (41) to (43) is solved to obtain the corresponding velocity pulse control quantity, so that the mission satellite is exactly at the orbital perigee position after the orbital plane transfer is completed. Δv pe is called the orbital-raising velocity increment. The mission satellite uses the difference between the perturbation forces received by it and the satellites in the target constellation to achieve transfer between adjacent orbital planes of the target constellation, without consuming any fuel for maneuvering, which significantly prolongs the on-orbit service life of the mission satellite.

[0139] The fsolve function in MATLAB is used to solve, and the maneuvering pulse required to be applied at time t2 is obtained.

[0140] Step three: Based on the tangential velocity pulse control quantity obtained in step two, a correction maneuver is applied during the orbital plane drift of the mission satellite to achieve matching of the perigee amplitude angle of the mission satellite and the phase of the target satellite. Repeat step one to calculate the relative position and velocity state of the mission satellite when it passes the target. Take the minimum correction maneuver velocity increment as the objective function, take the relative position of the target satellite that satisfies the minimum angular distance of the perigee when the mission satellite ends as the constraint, and take the time when the correction maneuver is applied as the optimization variable to construct an optimization problem. The relationship between the time when the correction maneuver is applied and the size of the correction maneuver velocity increment is established by the Lambert problem, and the optimization problem is solved to obtain the control quantity of the intermediate correction maneuver applied during the orbital plane transfer of the mission satellite. Since the mission satellite and several satellites on the adjacent orbital plane all have matching opportunities, under the premise of satisfying the constraint that the mission satellite is exactly at the perigee after the transfer is completed in step two, only the intermediate correction needs to be applied to achieve the coincidence of the mission satellite and a certain target satellite after the orbital plane transfer is completed.

[0141] The intermediate correction velocity increment Δv inter is applied in the manner of satisfying the following terminal constraint: at time t3, the orbital perigee of the mission satellite coincides with a certain target satellite. The terminal constraint is expressed as

[0142]

[0143] In the formula: ω c3 represents the perigee amplitude angle of the mission satellite at time t3, and ui denotes the latitude amplitude of the i-th target star on the 2nd orbit plane, and N denotes the number of target stars on the target orbit.

[0144] find the target star i with the minimum perigee distance to the mission star at t3 * , i.e.

[0145]

[0146] repeat step one to get the target orbit of the mission star with the i * -th target star as the starting point, which is called the target mission orbit. The perigee of the target mission orbit is the target position of the mission star at t3, denoted as

[0147] Consider placing the midcourse correction time t2 within the orbit of the mission star before reaching the target mission orbit. Partially optimize the correction impulse Δv inter to minimize the impact of the midcourse correction on the drift time of the orbit plane, i.e.

[0148] min J' = || Δv inter || (46)

[0149] satisfy the following constraints

[0150]

[0151] t3-t2≤T c1 (48)

[0152] where: r c3 denotes the position of the target star at t3, T c1 denotes the orbital period of the mission star after applying Δv pe . The Gauss algorithm is used to solve the Lambert problem, and the fmincon function can be used to solve the above partial optimization problem.

[0153] Step four: when the mission star reaches the relative position of the target star obtained in step three, apply a pulse Δv pe2 to make the relative velocity of the mission star satisfy the flyby condition, so that the mission star enters the resonance orbit of the flyby of the satellite in the star constellation, and the applied velocity pulse is recorded as the on-orbit velocity increment control quantity.

[0154] Step five: generate the control timing instructions of the task star to traverse all the satellites of the adjacent orbital planes of the target constellation in sequence according to the time sequence of the ascending track speed increment obtained in step two, the modified speed increment obtained in step three and the entry track speed increment obtained in step four, so as to realize the spacecraft trajectory planning of the constellation multi-target traversal flyby. According to the trajectory planning result of the task star, the fuel consumption required by the task star to traverse the targets can be reduced, and additional waiting adjustment phase is not required, so that the execution efficiency of the patrol task is improved.

[0155] To verify the feasibility of the method, it is assumed that there is a target constellation with Walker-delta as the basic configuration, and a spacecraft trajectory planning method for constellation multi-target traversal flyby is used to realize flyby traversal of all target satellites on the adjacent two orbital planes.

[0156] According to the above steps, the time sequence of the task star executing the maneuver is planned, a total of 3 maneuvers are executed, respectively after flying by the last target satellite of the first orbital plane, within the first circle of the orbital plane drift completion, and at the perigee position after the orbital plane drift completion, and the corresponding speed increment sizes are 141.1 m / s, 80.8 m / s and 149.8 m / s respectively, and the total speed increment consumed is 371.7 m / s. The drift time of the task star between the adjacent orbital planes is 136.2 hours, which is about 5.7 days, and the time for completely traversing the target stars on each orbital plane is 36.67 hours, so a total of 8.7 days is required to complete the traversal flyby of all targets on the two orbital planes. The target star RTN coordinate system established in this embodiment is shown in Figure 2 . The relative motion trajectory of the task star when flying by the target star is shown in Figure 2 . The coordinate system of the relative motion trajectory is the RTN coordinate system in Figure 3 . The green curve in the figure is the relative motion trajectory of the task star, the red dot is the mass center position of the target spacecraft, and the colored ball represents the flyby relative distance constraint condition. It can be seen that the green trajectory curve is tangent to the ball, and the constraint condition is satisfied. The state of the task star when transferring from the first target orbital plane to the second target orbital plane is shown in Figure 4 . The blue dashed line represents the initial trajectory of the target star before the orbital plane drift, and the blue solid line represents the task star orbit after the J2 perturbation drift. The task star is represented by a blue dot, and the target star is represented by a red dot. It can be seen that at this moment, the task star almost coincides with one of the target stars, and enters the resonance orbit with a semi-major axis slightly longer than the target star, indicating that the task star transfer rendezvous across the orbital plane achieves the expected design purpose.

[0157] The above detailed description of the specific description, the purpose, technical scheme and beneficial effects of the application are further described in detail, it should be understood that the above description is only a specific embodiment of the present application, for explaining the present application, and is not used to limit the protection scope of the present application, any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A method for spacecraft trajectory planning for constellation-oriented multi-target flyby, characterized in that: Comprising the following steps, Step one: taking the relative position and velocity state of the mission star when flying by the target star as the optimization variable, introducing the flyby time constraint and the ideal flyby azimuth, reducing the number of actual optimization variables by using the geometric relationship constraint of the flyby time and the orbital resonance relationship between the mission star and the target star, making the established mission star flyby state optimization problem contain only 3 degrees of freedom, converting the optimization of the state vector into the optimization of the space angle; relaxing the ideal flyby azimuth constraint as the minimum objective function, establishing the space angle value range as the constraint condition, constructing the mission star flyby state optimization problem, solving the relative position and velocity state to generate the relative motion trajectory of the mission star flying by a single target; due to the introduction of the resonance orbit constraint condition, the target star can naturally realize the traversal flyby of all target stars on the same orbital plane on the orbit of flying by a single target star; Step two: taking the time when the mission star flies by the last target on the single orbital plane of the star cluster according to the relative motion trajectory optimized in step one as the transfer departure time, adjusting the semi-major axis by applying a tangential velocity pulse, and realizing the orbital plane transfer of the mission star by using the difference in the right ascension drift rate of the mission star and the target star under J2 perturbation at different orbital semi-major axes; in order to realize that the mission star is exactly at the perigee position after completing the orbital plane transfer, according to the relationship that the drift transfer time is an integer multiple of the mission star period, an equation is constructed with the tangential velocity pulse applied at the transfer departure time as the only independent variable, and the corresponding velocity pulse control quantity is solved to realize that the mission star is exactly at the perigee position after completing the orbital plane transfer; the mission star realizes the transfer between adjacent orbital planes in the target star cluster by using the difference in the perturbation forces it and the satellites in the target star cluster receive; Step three: based on the tangential velocity pulse control quantity at the transfer departure time obtained in step two, a correction maneuver is applied during the orbital plane drift of the mission star to realize the matching of the mission star's perigee amplitude angle and the target star's phase; step one is repeated to calculate the relative position and velocity state of the mission star when flying by the target; taking the minimum correction maneuver velocity increment as the objective function, taking the relative position of the target star that satisfies the minimum flyby perigee angular distance of the mission star's terminal position as the constraint, and taking the time when the correction maneuver is applied as the optimization variable, an optimization problem is constructed; the relationship between the time when the correction maneuver is applied and the size of the correction maneuver velocity increment is established by using the Lambert problem, and the optimization problem is solved to obtain the control quantity of the intermediate correction maneuver applied during the orbital plane transfer of the mission star; since the mission star and several satellites on the adjacent orbital planes all have matching opportunities, only the intermediate correction needs to be applied to realize that the mission star coincides with a certain target star after completing the orbital plane transfer under the premise that the mission star is exactly at the perigee after the transfer is completed in step two; Step four: when the mission star reaches the relative position of the flyby target star obtained in step three, apply a pulse Δv pe2 The relative velocity of the mission star satisfies the traversal flyby condition, the mission star enters the resonance orbit of the satellite traversal flyby in the star cluster, and the applied velocity pulse is recorded as the on-orbit velocity increment control quantity; Step five: sequentially connect the ascending track speed increment obtained in step two, the modified speed increment obtained in step three and the entry track speed increment obtained in step four in time sequence to generate control time sequence instructions for the task star to traverse all the satellites on the adjacent orbital planes of the target constellation, so as to realize the spacecraft trajectory planning for the constellation multi-target traversal; and the guidance control of the task star according to the trajectory planning result of the task star can reduce the fuel consumption required by the task star for traversing and patrolling between targets, and does not need to perform additional waiting adjustment, thereby improving the execution efficiency of the patrol task.

2. The method of claim 1, wherein: Step one is implemented by, Step 1.1 establishes a target star RTN coordinate system, and establishes a general expression of the position of the target star in terms of the classical orbital elements according to the characteristics of the multi-target star co-orbital plane of the Walker-δ constellation configuration; An RTN coordinate system is established with the target star as the origin, the x-axis points to the outside along the radial direction of the target star, the y-axis is perpendicular to the x-axis in the orbital plane of the target star and points to the velocity direction, and the z-axis is perpendicular to the orbital plane of the target star and forms a right-hand system with the x-axis and the y-axis; In the Walker-δ constellation configuration, the orbital elements of the target stars in the same orbital plane only differ in the true anomaly, and the orbital elements of the jth target satellite are represented by a vector as shown in formula (1) where: a t denotes the semi-major axis of the orbit, e t denotes the eccentricity, i t denotes the inclination of the orbit, Ω t denotes the longitude of the ascending node, ω t denotes the argument of the perigee, θ0denotes the true anomaly at the initial time, N denotes the total number of satellites in this orbit, for Step 1.2 Take the time when the distance between the target star and the mission star is the shortest as the flyby time t s , and take the relative position vector δr s and the relative velocity vector δv s of the mission star and the target star at this time as the optimization variables; According to the load capacity of the task star and the mission purpose, the relative position and relative velocity constraints of the task star when flying by the target are set; Setting the flyby time task star state δr s and δv s satisfy the following constraint conditions δr s = ||δr s || = r fix (2) δv s =||δv s ||∈[v L ,v U ) (3) where: r fix is the fixed distance of the task star flying over the target star, v L , v U are the lower bound and upper bound of the flying-over relative velocity, respectively. The ideal orientation of the mission star at the time of the flyby of the target is indicated by the azimuth angle α t and the elevation angle β t ; the corresponding unit vectors are denoted by l tar ; The ideal flyby position vector [l tar ] t As follows [l tar ] t = [cosβ t cosα t cosβ t sina t sinβ t ] T (4) Step 1.3 utilizes the constraints between the mission satellite and the target at the flyby moment to reduce and transform the optimization variables described in Step 1.2, and then transforms δr. s As a spatial vector, δv is described by azimuth angle α and elevation angle β. s Then, as perpendicular to δr s A vector in a plane is represented by its direction angle γ and relative velocity magnitude δv. s The description proceeds as follows: the resonance ratio between the mission satellite's orbit and the target satellite's orbit is determined by the number of satellites contained in each orbital plane of the target constellation, thereby determining the semi-major axis of the mission satellite's orbit. The vitality equation is then used to construct its relationship with the relative velocity δv. s The relationship between these factors further reduces the number of actual optimization variables, resulting in a mission star flyby state optimization problem with only 3 degrees of freedom. The ideal flyby orientation constraint established in step 1.2 is relaxed into an objective function, with the actual flyby position δr as the objective function. s and ideal flight direction l tar The goal is to minimize the included angle between them, thus completing the construction of the task star's flight state optimization problem. The position and velocity state vector of the target star in the inertial frame at time t is denoted as r s and v ts ; the relative position [δr ts ] s is expressed in the target star RTN coordinate system using two angles α and β, and the relative velocity [δv s ] t is expressed as [δr s ] t = δr s • [cos β cos α cos β sin α sin β] T (5) At the same time, the minimum distance necessary condition as shown in formula (6) is satisfied to ensure that the task star and the target star are at the closest distance δr s • δv s = 0 (6) δv = δv + δv s Applying additional constraints, its coordinate components in the RTN coordinate system [δv s ] t is expressed as [δv s ] t = δv s L z (-α)L y (β)[0 cosγ sinγ] T (7) where L x , L y , L z are rotation matrices along the x, y, z axes; Tracking stars in t s The state at time should be wherein r cs , r ts represent the positions of the mission star and the target star at time t s , respectively, v cs , v ts represent the velocities of the mission star and the target star, respectively, [δr s ] eci , [δv s ] eci are the components of the relative position and relative velocity vectors in the Earth-Centered Inertial (ECI) frame, obtained by coordinate transformation from equations (5) and (7) [δr s ] eci = L z (-Ω t )L x (-i t )L z (-θ t -ω t )[δr s ] t (9) [δv s ] eci = L z (-Ω t )L x (-i t )L z (-θ t -ω t )[δv s ] t (10) To meet the condition that the mission star still passes the (j+1)th target star in the same relative state, the orbit period T of the mission star should satisfy c should satisfy where T represents the target satellite orbit period, N represents the total number of satellites in the orbit plane, and a represents the semi-major axis of the tracking satellite calculated according to the relationship between the semi-major axis and the orbit period t where T represents the target satellite orbit period, N represents the total number of satellites in the orbit plane, and a represents the semi-major axis of the tracking satellite calculated according to the relationship between the semi-major axis and the orbit period c is In formula, μ is the gravitational constant of the earth; according to the energy equation, formula (8) is satisfied wherein: r cs = ||r cs ||, v cs = ||v cs ||; The performance index function is defined as the angle between the ideal direction vector l tar and the actual flyover bearing [δr s ] t / ||δr s ||, denoted as J = min [arccos(l tar • δr s / ||δr s ||)] (14) According to the above analysis and modeling process, the state optimization problem of the task star flying by the target star has the following four free variables One constraint formula (13) is satisfied to constitute a 3-degree-of-freedom task star flying state optimization problem; Step 1.4 solves the task star flying state optimization problem established in step 1.3 to obtain the relative state of the task star when flying by the target star, and constitutes the trajectory of the task star when flying by a single target star; since the resonance orbit constraint condition is introduced, the target star can naturally realize the traversal flyby of all the target stars on the same orbital plane on the orbit of a single target star.

3. The method of claim 2, wherein: Step two is implemented by, The time when the task star applies an acceleration pulse at the perigee is recorded as t2, and the time when the task star reaches the adjacent orbital plane and flies by the first target star is recorded as t3; At t3, the task star should be at the perigee of the orbit, that is, the orbital plane drift time of the task star is an integer multiple of the task star period; the above integer multiple relationship is represented as where mod(·) function denotes the remainder, ΔT is the orbital plane transfer time obtained in equation (19), T c1 denotes the orbital period of the mission satellite after the pulse is applied at time t1; where ΔT is given by the difference between the ascending node rates of the mission satellite and the target satellite, under the influence of J2 perturbation, the ascending node right ascension will drift according to the Gauss perturbation equation as shown in equation (17), to obtain the ascending node right ascension rate of change of the mission satellite is: where: a c denotes the semi-major axis of the target star, i c denotes the inclination of the target star orbit; the rate of change of the right ascension of the ascending node of the target star can be obtained in the same way The expression for From equation (17) and (18), the orbital plane of the mission satellite can be changed by changing the semi-major axis a of the mission satellite c The orbital plane of the mission satellite can be changed by changing the difference between the ascending node velocities of the mission satellite and the target satellite. The difference in right ascension of the ascending nodes between adjacent orbital planes of the target satellite is denoted as ΔΩ. The time required to complete the drift is The tangential velocity increment applied at time t2is denoted as Δv pe Substituting equations (17) to (19) into equation (16) gives an equation with Δv pe as the only unknown. Solving the equation group composed of equations (16), (17) to (19) gives the corresponding velocity impulse control quantity, which makes the mission satellite just at the orbit perigee position after the orbit plane transfer is completed, and Δv pe is called the orbit-raising velocity increment; the mission satellite uses the difference between the perturbation forces it and the satellites in the target constellation receive to realize the transfer between adjacent orbit planes in the target constellation.

4. The method of claim 3, wherein: Step three is implemented by, The terminal constraint that the mission satellite orbit apsis coincides with a target satellite at time t3 is satisfied by applying a midcourse correction velocity increment Δv inter The terminal constraint is expressed as ω c3 = u i , wherein ω c3 denotes the argument of perigee of the mission satellite at time t3, u i denotes the latitude argument of the i-th target satellite on the 2nd orbital plane, N denotes the number of target satellites on the target orbit; Finding the target star i with the minimum perigee distance angle from the task star at time t3 * i.e. Step one is repeated to obtain a flyby mission orbit with the target star numbered i * as the starting point, as the target mission orbit of the target star. The near apsis of the target mission orbit is the target position of the mission satellite at time t3, denoted as Considering that the midcourse correction time t2 is placed within one revolution of the mission star to the target mission orbit, the correction impulse Δv inter Local optimization is performed to minimize the impact of the midcourse correction on the drift time of the orbital plane, i.e. min J' = ||Δv inter || (22) The following constraints are satisfied t3 - t2≤ T c1 (24) where: r c3 represents the target star position at t3, T c1 represents the orbit period after the mission star applies Δv pe The Gauss algorithm is used to solve the Lambert problem, and the fmincon function is combined to solve the above local optimization problem.

5. The method according to claim 3 or 4, wherein: The fmincon function in MATLAB is used to solve the above optimization problem; The fsolve function in MATLAB is used to solve, and the maneuver pulse required at t2 is obtained.

Images