CN103488830B - The task simulation system that a kind of ground based on Cycler track moon comes and goes - Google Patents

Abstract

The invention discloses the task simulation system come and gone by a kind of ground based on Cycler track moon, coming and going between the earth and the moon of task can be emulated, across spacecraft orbit design, spacecraft task analysis and spacecraft dynamics emulation etc. by this system.The present invention comes resonance type Cycler track and chummage type Cycle Track desigh correction, Lambert transfer orbit design modification, Launch Encounter window analysis and the technical problem of ground moon shuttle system task analysis by building resonance Cycler model trajectory, structure chummage Cycler model trajectory, multiple shooting method orbital exponent module and Lambert transfer orbit acquisition module.Present system considers solar gravitation perturbation when obtaining ground moon Cycler track, and uses multiple shooting method to obtain approximate period Cycler track, and precision increases compared with the cycle Cycler track that existing method calculates.

Description

Earth-moon round-trip task simulation system based on Cycler orbit

Technical Field

The invention relates to a round-trip task simulation between the earth and the moon, in particular to a round-trip task simulation system based on the cycle orbit.

Background

The first lunar exploration project in China consists of five systems of a lunar exploration satellite, a carrier rocket, a launching field, measurement and control, ground application and the like, and the lunar exploration satellite project in China also has five project targets: firstly, the first lunar exploration satellite in China is developed and launched; secondly, the basic technology of lunar exploration is mastered preliminarily; thirdly, carrying out lunar scientific exploration for the first time; fourthly, primarily constructing a lunar exploration aerospace engineering system; fifthly, experience is accumulated for subsequent engineering of lunar exploration. Therefore, the key technology of lunar exploration satellites is broken through; preliminarily establishing a large system of Chinese deep space exploration engineering; verifying each key technology such as effective load and data interpretation; preliminarily establishing a Chinese deep space exploration technology development system; and culturing the corresponding talent team.

The first four scientific tasks of the lunar exploration project in China are as follows:

the method comprises the steps of acquiring a three-dimensional stereoscopic image of the lunar surface, finely dividing a basic structure and a landform unit of the lunar surface, researching the form, the size, the distribution, the density and the like of a collision pit on the lunar surface, providing basic data for the division of the surface age of the geostationary and the early evolution historical research, providing basic data for the location selection of a lunar surface soft landing area and the position optimization of a lunar base and the like.

And secondly, analyzing the distribution characteristics of the useful element content and the material type on the lunar surface, mainly surveying the content and the distribution of 14 elements such as titanium, iron and the like with development and utilization values on the lunar surface, drawing a full lunar distribution diagram of each element, a lunar rock, mineral and geological thematic map and the like, finding out the enrichment area of each element on the lunar surface, evaluating the development and utilization prospect of lunar mineral resources and the like.

And thirdly, detecting the thickness of lunar soil, namely acquiring the thickness data of lunar soil on the surface of the moon by using a microwave radiation technology so as to obtain the age and the distribution of the lunar surface, and estimating the content, the resource distribution, the resource quantity and the like of the nuclear fusion power generation fuel helium 3 on the basis.

And fourthly, detecting the space environment from the earth to the moon. The average distance between the moon and the earth is 38 kilometres, the satellite is located in a far magnetic tail region of the earth magnetic field space, the satellite can detect solar cosmic ray high-energy particles and solar wind plasmas in the region, and interaction between solar wind and the moon and interaction between the earth magnetic field magnetic tail and the moon are researched.

The Cycler orbit is an orbit that periodically travels between the earth and the moon, and flies around the planet without staying. An aircraft operating on a cyclic orbit can maintain inter-planet flight for long periods of time (years or even decades) without orbital maneuver (or with only minimal orbital maneuver), and is therefore considered an economical long-term mission for energy conservation. The Cycler orbit scheme can be further subdivided into a whole circle orbit, a half circle orbit, a common orbit and the like according to different numbers of turns around a day. The Cycler orbit can be divided into a resonance type Cycler orbit and a homosomic type Cycler orbit.

The Circular Restricted Three-Body Problem (CR 3BP) describes the movement of a third Body of relatively infinite mass under the influence of the attractive force of two main celestial bodies moving in a circle around their common centroid. Refer to the related contents on pages 19 to 21 of "doctor's academic thesis of research institute of Chinese academy of sciences" of Li Ming Tao, 5.2010.

The Bi-circular Model (BCM) is a basic Model for studying the motion law of an infinite small mass body under the action of the universal gravitation of a system of moon (running on a circular orbit around the system centroid), sun and earth (running on a circular orbit around a common centroid).

Disclosure of Invention

The invention aims to provide a earth-moon round-trip task simulation system based on a Cycler orbit, which generates a resonance Cycler orbit and a homoclinic Cycler orbit under a circular restrictive three-body model. Due to the eccentricity of the lunar orbit and the influence of the gravity of other celestial bodies (such as the sun, the wooden star and the like), the orbits established under the circular restrictive three-body model are different from the real situation, and even due to the influence of perturbation, the orbits cannot be kept constant. Therefore, the orbit established under the circular restrictive trisomy model is used as an initial value, optimization is carried out under the double-circle model, and the resonance type periodic orbit and the homoclinic type periodic orbit are corrected by using a multiple targeting method. The modified resonance type periodic orbit and the homoclinic type periodic orbit can meet the earth and the moon regularly, and a transmitting window which takes the earth or the moon as a target is generated. In order to complete the complete earth-moon transfer orbit, the simulation system of the invention also generates an 'earth-cycle transfer orbit' for earth takeoff and moon landing according to a classical solution of the two-body Lambert problem; "cycle transfer orbit-moon" for moon takeoff and earth landing; and modified under a four-body model. The fuel consumption and the flight time of the two transfer orbits were calculated and compared in the simulation system of the present invention.

The invention relates to a earth-moon round-trip task simulation system based on a Cycler orbit, which comprises a resonance Cycler orbit model 10, a co-homed Cycler orbit model 20, a multi-targeting method orbit correction module 30 and a Lambert transfer orbit acquisition module 50.

The first aspect of the method for constructing the resonant Cycler orbit model 10 is to construct M according to a circular restrictive three-body problem model CR3BP_CR3BPA kinetic model; in the second aspect, the double-circle model BCM is adopted to match the M_CR3BPOptimizing the dynamic model to obtain M_BCMAnd (4) a dynamic model.

The construction of the co-homed Cycler orbit model 20 is based on the construction of the M of the double-circle model BCM_LAnd (4) a dynamic model.

The multiple targeting trajectory modification module 30, in a first aspect, uses multiple targeting to respectively target M_BCMCorrecting to obtain a corrected resonance Cycler orbit dynamics model DM_BCM(ii) a Second aspect to DM_BCMMethod for obtaining resonant Cycler orbit SM by adopting fourth-order Runge Kutta method_BCM(ii) a In the third aspect, multiple targeting methods are used to target M separately_LCorrecting to obtain a corrected compatible Cycler orbit dynamics model DM_L(ii) a Fourth aspect to DM_LObtaining co-sited Cycler orbit SM by adopting four-step Runge Kutta method_L。

The Lambert transition track obtaining module 50 adopts a gaussian-global variable composite algorithm to pair SM in the first aspect_BCMProcessing to obtain a first set of Lambert transfer tracksThe second aspect adopts a differential correction algorithm pairCorrecting to obtain a first set of corrected Lambert transfer tracksThird aspect is to measure SM according to variation rule of emission point latitude and white-red intersection angle_BCMAnalyzing a launching window and an intersection window to obtain the launching and orbit entering time of the spacecraft; the fourth aspect adopts a Gaussian-global variable composite algorithm to SM_LProcessing to obtain a second set of Lambert transfer tracksThe fifth aspect adopts a differential correction algorithm pairCorrecting to obtain a second set of corrected Lambert transfer tracksThe sixth aspect can learn the SM by comparing the fuel consumption amount and the total flight time of the earth-moon round trip mission_BCMAnd SM_LThe respective advantages of the orbits provide optimized design indexes for the earth-moon round trip task of the deep space exploration of the low-thrust spacecraft.

The simulation system of the invention has the advantages that:

firstly, the invention generates a resonance type Cycler orbit and a homoclinic type Cycler orbit under a circular restrictive trisomy model. The method aims to correct the influence of solar attraction perturbation on track design in the existing resonant type Cycler track and the homoclinic type Cycler track.

Secondly, the system considers the solar attraction perturbation when obtaining the earth-moon Cycler orbit, obtains the approximate cycle Cycler orbit by using a multiple targeting method, and improves the precision compared with the cycle Cycler orbit calculated by the existing method.

The system provides an iterative algorithm when analyzing the intersection window of the spacecraft from the ground to the circulator orbit and the intersection window of the spacecraft from the moon to the resonance circulator orbit, can obtain the transfer window under the condition of considering the solar attraction perturbation, and has higher precision compared with the calculation result of the existing Lambert method.

Drawings

FIG. 1 is a block diagram of the construction of a earth-moon round trip task simulation system based on the Cycler orbit.

FIG. 2 is a schematic diagram of the coordinate system of Earth-moon centroid inertial coordinate system O-XYZ and Earth-moon centroid rotational coordinate system O-XYZ.

FIG. 2A is a day-to-earth rotation coordinate system O_S-X_SY_SZ_SSchematic diagram of the coordinate system of (1).

FIG. 3 is a diagram of a resonance Cycler orbital dynamics model DM modified by the multiple targeting method of the present invention_BCMFigure (a).

FIG. 3A is a modified homeostatic Cycler orbital dynamics model DM by the multiple targeting method of the present invention_LFigure (a).

Fig. 4 is a schematic diagram of two-point boundary value acquisition in space dynamics under the Lambert problem.

FIG. 5 is a diagram of SM processed by the Gaussian-global variable composite algorithm of the present invention_BCMAnd (5) drawing a simulation result.

FIG. 5A is a graph of simulation results for a first set of modified Lambert transition tracks after modification by a differential correction algorithm.

FIG. 5B is a graph of simulation results for a second set of modified Lambert transition tracks after modification by the differential correction algorithm.

Fig. 6 is a schematic view of the window for direct track entry of the west chang.

FIG. 6A is a schematic of the trajectory of the velocity increment required for the most fuel efficient Hotman transfer.

FIG. 6B is a schematic of the trajectory of the velocity increment required for the most fuel Hotman transfer.

FIG. 6C is a schematic diagram of the emission trajectory under Earth-moon centroid inertial coordinate system O-XYZ.

FIG. 6D is a schematic diagram of the emission trajectory under the Earth-moon centroid rotating coordinate system O-xyz.

Fig. 7 is a schematic view of the window for weschang direct orbit entry for lunar takeoff.

Fig. 7A is a schematic diagram of the search result under the constraint that the height of the parking track is 100 km.

Fig. 7B is a schematic diagram of an entry window having a geographic longitude of [ -170 °, -135 ° ] [ -43 °, -5 ° ] [ -45 °,108 ° ] [150 °,165 ° ].

Fig. 7C is a graph of results of simulation of the co-localized Lambert tracks for the two-body Lambert problem.

FIG. 7D is a second set of modified Lambert transition tracks after a differential correction algorithmAnd (4) obtaining a simulation result.

Detailed Description

The present invention will be described in further detail with reference to the accompanying drawings.

Referring to fig. 1, the invention is a earth-moon round-trip task simulation system based on a Cycler orbit, which comprises a resonant Cycler orbit model 10, a co-homed Cycler orbit model 20, a multi-targeting method orbit correction module 30 and a Lambert transfer orbit acquisition module 50.

The first aspect of constructing the resonant Cycler orbit model 10 is based on circular constraintsConstruction of M by sex trisomy problem model CR3BP_CR3BPA kinetic model; in the second aspect, the double-circle model BCM is adopted to match the M_CR3BPOptimizing the dynamic model to obtain M_BCMAnd (4) a dynamic model.

The construction of the co-homed Cycler orbit model 20 is based on the construction of the M of the double-circle model BCM_LAnd (4) a dynamic model.

The multiple targeting trajectory modification module 30, in a first aspect, uses multiple targeting to respectively target M_BCMCorrecting to obtain a corrected resonance Cycler orbit dynamics model DM_BCM(ii) a Second aspect to DM_BCMMethod for obtaining resonant Cycler orbit SM by adopting fourth-order Runge Kutta method_BCM(ii) a In the third aspect, multiple targeting methods are used to target M separately_LCorrecting to obtain a corrected compatible Cycler orbit dynamics model DM_L(ii) a Fourth aspect to DM_LObtaining co-sited Cycler orbit SM by adopting four-step Runge Kutta method_L。

The Lambert transition track obtaining module 50 adopts a gaussian-global variable composite algorithm to pair SM in the first aspect_BCMProcessing to obtain a first set of Lambert transfer tracksThe second aspect adopts a differential correction algorithm pairCorrecting to obtain a first set of corrected Lambert transfer tracksThird aspect is to measure SM according to variation rule of emission point latitude and white-red intersection angle_BCMAnalyzing a launching window and an intersection window to obtain the launching and orbit entering time of the spacecraft; the fourth aspect adopts a Gaussian-global variable composite algorithm to SM_LProcessing to obtain a second set of Lambert transfer tracksThe fifth aspect adopts a differential correction algorithm pairCorrecting to obtain a second set of corrected Lambert transfer tracksThe sixth aspect can learn the SM by comparing the fuel consumption amount and the total flight time of the earth-moon round trip mission_BCMAnd SM_LThe respective advantages of the orbits provide optimized design indexes for the earth-moon round trip task of the deep space exploration of the low-thrust spacecraft.

(I) construction of a resonant Cycler orbit model 10

In the invention, a resonance Cycler orbit dynamics model M is obtained according to a circular restrictive three-body problem model CR3BP_CR3BP：

$M_{CR3BP}=\left\{\begin{array}{c}\frac{\∂U}{\∂x}={\mathrm{HA}}^x-2{\mathrm{VA}}^y\\ \frac{\∂U}{\∂y}={\mathrm{HA}}^y+2{\mathrm{VA}}^x\\ \frac{\∂U}{\∂z}={\mathrm{HA}}^z\end{array}\right.---\left(1\right)$

As shown in FIG. 2, Earth-moon centroid inertial coordinate system O-XYZ and Earth-moon centroid rotational coordinate system O-XYZ. Wherein O is the earth-moon mass center, and the earth-moon system motion plane is an XY coordinate plane. In the earth-moon mass center inertial coordinate system, X is the direction in which the earth points to the moon at the initial moment, Z is the angular velocity direction of the earth-moon system, and Y is determined according to the X and Z right-hand rules; the Z direction in the earth-moon mass center rotating coordinate system is coincided with the Z direction of the earth-moon mass center inertial coordinate system, the x direction is always the direction in which the earth points to the moon, and the y is established according to the right-hand system rule. Let the earth be P₁The moon is P₂The spacecraft is P, and the position of the spacecraft P under the earth-moon mass center rotating coordinate system O-xyz is recorded as (x)_P,y_P,z_P) (ii) a Spacecraft P to earth P₁Is recorded asSpacecraft P to moon P₂Is recorded asThe distance from the spacecraft P to the Earth-moon centroid O is denoted as R.

Under the earth-moon mass center rotating coordinate system O-xyz, the physical meaning of the letter of formula (1) is as follows:represents the partial derivative in the x-direction;

u represents the potential function of the spacecraft P, andμ₁represents the mass ratio of the moon to the earth, generally taking the value of 0.01215;

represents the partial derivative in the y-direction;

represents the partial derivative in the z direction;

VA^xrepresents the velocity of the spacecraft P in the x direction;

VA^yrepresents the velocity of the spacecraft P in the y direction;

HA^xrepresents the acceleration of the spacecraft P in the x direction;

HA^yrepresents the acceleration of the spacecraft P in the y direction;

HA^zrepresenting the acceleration of the spacecraft P in the z direction.

In the present invention, according to the double round model BCM, M is paired_CR3BPCorrecting to obtain a dynamic model M under the double-circle model BCM_BCM：

$M_{BCM}=\left\{\begin{array}{c}{\mathrm{HB}}^{X_S}=2{\mathrm{VB}}^{Y_S}+{\mathrm{xs}}_P-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{xs}}_P+{\mathrm{\μ}}_2}{r_{PS}^3}-\\ {\mathrm{\μ}}_2\frac{{\mathrm{xs}}_P-1+{\mathrm{\μ}}_2}{r_{PE}^3}-m_M\frac{{\mathrm{xs}}_P-x_{P_2}}{r_{PM}^3}\\ {\mathrm{HB}}^{Y_S}=-2{\mathrm{VB}}^{X_S}+{\mathrm{ys}}_P-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{ys}}_P}{r_{PS}^3}-\\ {\mathrm{\μ}}_2\frac{{\mathrm{ys}}_P}{r_{PE}^3}-m_M\frac{{\mathrm{ys}}_P-y_{P_2}}{r_{PM}^3}\\ {\mathrm{HB}}^{Z_S}=-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{xs}}_P}{r_{PS}^3}-{\mathrm{\μ}}_2\frac{{\mathrm{zs}}_P}{r_{PE}^3}-m_M\frac{{\mathrm{zs}}_P}{r_{PM}^3}\end{array}\right.---\left(2\right)$

In the double circle model BCM, a day-ground rotation coordinate system O is used_S-X_SY_SZ_SAs shown in FIG. 2A, the sun is P₃The spacecraft P rotates the coordinate system O in the sun and the earth_S-X_SY_SZ_SThe lower position is denoted as (xs)_P,ys_P,zs_P) (ii) a Moon P₂In the sun-ground rotation coordinate system O_S-X_SY_SZ_SThe lower position is notedr_PSRepresents the dimensionless distance of the spacecraft P from the sun, anr_PERepresents the dimensionless distance of the spacecraft P from the earth, anr_PMRepresents the dimensionless distance of the spacecraft P from the moon, anμ₂Representing the mass ratio of the earth to the sun, typically 0.000003003, m_MRepresenting the dimensionless quality of the moon. By the centroid of the sun and earth_SIs an origin, and the sun-ground connecting line is X_SAxis, pointing towards the earth, X_SPositive direction of the axis. Y is_SThe axis being perpendicular to the plane of motion of the sun, Z being determined according to the right-hand rule_SThe direction of the axis.

In the sun-ground rotation coordinate system O_S-X_SY_SZ_SThe letter in the formula (2)The principle meaning is:

indicates M is_BCMSpacecraft under model P at X_SA speed in a direction;

indicates M is_BCMSpacecraft under model P in Y_SA speed in a direction;

indicates M is_BCMSpacecraft under model P at X_SAcceleration in a direction;

indicates M is_BCMSpacecraft under model P in Y_SAcceleration in a direction;

indicates M is_BCMSpacecraft under model P in Z_SAcceleration in the direction.

(II) construction of a co-sited Cycler orbit model 20

In the invention, a co-homed Cycler orbit dynamics model M is obtained according to a double-circle model BCM_L：

$M_L=\left\{\begin{array}{c}{\mathrm{HC}}^{X_S}=2{\mathrm{VC}}^{Y_S}+{\mathrm{xs}}_P-(1-\mathrm{\μ})\frac{{\mathrm{xs}}_P+\mathrm{\μ}}{r_{PS}^3}-\\ {\mathrm{\μ}}_2\frac{{\mathrm{xs}}_P-1+{\mathrm{\μ}}_2}{r_{PE}^3}-m_M\frac{{\mathrm{xs}}_P-x_{P_2}}{r_{PM}^3}\\ {\mathrm{HC}}^{Y_S}=-2{\mathrm{VC}}^{X_S}+{\mathrm{ys}}_P-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{ys}}_P}{r_{PS}^3}-\\ {\mathrm{\μ}}_2\frac{{\mathrm{ys}}_P}{r_{PE}^3}-m_M\frac{{\mathrm{ys}}_P-y_{P_2}}{r_{PM}^3}\\ {\mathrm{HC}}^{Z_S}=-(1-{\mathrm{\μ}}_2)\frac{{\mathrm{zs}}_P}{r_{PS}^3}-{\mathrm{\μ}}_2\frac{{\mathrm{zs}}_P}{r_{PE}^3}-m_M\frac{{\mathrm{zs}}_P}{r_{PM}^3}\end{array}\right.---\left(3\right)$

In the sun-ground rotation coordinate system O_S-X_SY_SZ_SThe physical meaning of the letters in formula (3) is:

indicates M is_LSpacecraft under model P at X_SA speed in a direction;

indicates M is_LSpacecraft under model P in Y_SA speed in a direction;

indicates M is_LSpacecraft under model P at X_SAcceleration in a direction;

indicates M is_LSpacecraft under model P in Y_SAcceleration in a direction;

indicates M is_LSpacecraft under model P in Z_SAcceleration in the direction.

(III) multiple targeting trajectory correction Module 30

In the present invention, multiple targeting is used to pair initialized M_BCMThe model is corrected to obtain a corrected resonance Cycler orbit dynamics model DM_BCM. The DM_BCMCharacterized in a graph (as shown in fig. 3). The multiple targeting method is described in "calculation of optimal flight trajectory for a class", volume ninth, first phase, pages A21-A22, written in 1988, by Wangbender et al.

M in one cycle_BCMTaking 500 samples at equal time intervals on the model, and marking any sample as ini _ b (i) ((i))1, …, 500); according to M_BCMThe model is integrated by a fourth-order Runge Kutta method, and the position and speed state quantity of any point is recorded as ini _ a (i) (1, … and 500), wherein ini _ a (i) is a matrix of 500 × 6.

At each time interval, 100 samples are taken and integrated with ODE45, taking the position velocity state quantity φ (ini _ a (i): at the end point). Note that the difference between the position-velocity state quantity phi (ini _ a (i): at the end point) and the position-velocity state quantity ini _ a (i:): at any point is F (i): phi (ini _ a (i):)) -ini _ a (i:). The difference F (i:) is applied to allow the calculated values and the actual calculated values to be error free, and should be minimized as much as possible during the simulation.

For the above difference F (i,: minimum value is obtained by using newton iteration method, which refers to MATLAB implementation of newton iteration method published in 2011, sixth stage, page 20, author yunlei, where the set precision is 1 × 10^-10Under the condition, the set precision is achieved through 8 iterations.

Iteration 1	2.4925020681966167×10-1
		Iteration 2	9.5043845337684230×10-2
Iteration 3	2.3381313323453772×10-2
		Iteration 4	1.9075592366893615×10-2
Iteration 5	4.5915032506800547×10-4
		Iteration 6	8.9114511750870658×10-5
Iteration 7	4.5782899616269825×10-9
		Iteration 8	8.9656802854146783×10-10

If the accuracy is larger than the set accuracy, repeating the iteration by using a Newton iteration method until the accuracy is smaller than the set accuracy. Thus, the corrected model DM can be obtained_BCM。

For DM_BCMMethod for obtaining resonant Cycler orbit SM by adopting fourth-order Runge Kutta method_BCMThe SM of_BCMThe track is shown in figure 3. In the figure, each point reflects the position of the spacecraft P under the earth-moon mass center inertial coordinate system O-XYZ in five cycles. The fourth-order Runge Kutta method refers to the Runge Kutta method and its Mathemica implementation published in 2006, Vol.18, No. 2, pages 72 to 73, author Chen 35468, and allergy.

In the present invention, multiple targeting is used to pair initialized M_LThe model is corrected to obtain a corrected compatible Cycler orbit dynamics model DM_L. The DM_LCharacterized in a graph (as shown in fig. 3A).

M in one cycle_LSampling 500 samples at equal time intervals on the model, and recording any sample as ini _ d (i) (1, … and 500); according to M_LThe model is integrated by adopting a four-step Runge Kutta method to obtain any pointThe position and speed state quantities are represented as ini _ c (i) (1, …,500), and ini _ c (i, 500) is a matrix of 500 × 6.

At each time interval, 100 samples are taken and integrated with ODE45, and the position velocity state quantity θ (ini _ c (i:)) of the end point is extracted. The difference between the position-speed state quantity θ (ini _ c (i): at the end point) and the position-speed state quantity ini _ c (i:): at any point is G (i): i.e.,: θ (ini _ c (i):)) -ini _ c (i:). The difference G (i:) is applied so that the calculated values and the actual calculated values are error free and should be minimized as much as possible during the simulation.

For the above difference G (i:) a Newton iteration method is used to find the minimum value with a set accuracy of 1 × 10^-10Under the condition, the set precision is achieved through 8 iterations.

If the accuracy is larger than the set accuracy, repeating the iteration by using a Newton iteration method until the accuracy is smaller than the set accuracy. Thus, the corrected model DM can be obtained_L。

For DM_LMethod for obtaining resonant Cycler orbit SM by adopting fourth-order Runge Kutta method_LThe SM of_LThe track is shown in fig. 3A. In the figure, each point reflects the position of the spacecraft P under the Earth-moon mass center inertial coordinate system O-XYZ.

For DM_LObtaining co-sited Cycler orbit SM by adopting four-step Runge Kutta method_L。

(IV) Lambert transition track acquisition module 50

The Lambert problem is a two-point boundary value problem in space dynamics, and is used for meeting of spacecrafts and interception of missilesAnd has wide application in the fields of interplanetary navigation and the like. As shown in fig. 4, the initial end point E of the spacecraft P₁Terminal E₂Respectively is L₁And L₂The focus of the ellipse transfer orbit is located at the geocentric P₁. Initial endpoint E₁Is recorded as t₁Terminal E₂Is recorded as t₂And θ is the transfer angle.

According to Lambert's law of flight time, the transfer time between any two points on the elliptical transfer orbit is equal to the sum (L) of the radius of the major semi-axis ra and the radius of the two points of the elliptical transfer orbit₁+L₂) And center angle θ, then expressed as:

t_f＝W(ra,(L₁+L₂),R_P) (4)

if L is₁And L₂The sum is constant, the major semi-axis ra is constant, and the initial endpoint E₁And terminal E₂A distance R therebetween_PIs constant, then from the initial end point E₁To terminal E₂Time of flight transition t_fIs also a constant. The determination of the elliptical transfer orbit and the choice of the speed between the two points are the key to the Lambert problem, which is described as a gaussian problem as follows: tracking spacecraft P departure position vector (L)₁) And velocity vector (v)₁) Spacecraft P terminal position vector (L)₂) And velocity vector (v)₂) Flight transition time of t_fThe spacecraft P in the elliptical transfer orbit E₁The initial velocity at a point is v₁₀The spacecraft P in the elliptical transfer orbit E₂The ending velocity at a point is v₂₀. The goal of the Gaussian problem is to solve for the initial position velocity delta Δ ν₁Velocity increment Δ v from terminal position₂And the magnitude of the applied pulse thrust is determined.

The gaussian problem can be solved by the following transcendental system of equations:

L₂＝k×L₁+g×v₁₀(5)

$v_{20}=\stackrel{\·}{k}\×L_1+\stackrel{\·}{g}\×v_{10}---\left(6\right)$

lagrange coefficient of one

Lagrange coefficient of two

Lagrange coefficient of three

Lagrange coefficient of fourMu is the mass ratio of the secondary celestial body to the main celestial body, and the value is 0.01215 in a circular restrictive three-body problem model CR3 BP; h is a variable on the elliptical transfer orbit of the spacecraft P at different positions.

In the formula (5) and the formula (6), L is known₁And v₁₀Or L is₂And v₂₀The elliptical transfer orbit followed by the spacecraft P can be determined. It is clear that once the Lagrangian coefficients k, g are determined,The gaussian problem can be readily solved.

In the invention, a global variable algorithm is adopted for solving the formula (5) and the formula (6).

The Gaussian-global variable composite algorithm refers to 'space rendezvous and docking mission planning' published in 2008, contents from pages 81 to 100, pages 142 to 144, authors in Tang national gold and the like.

Assuming a moon to resonant Cycler orbit transfer time of 3 days, the Lambert transfer orbit starting point was 100km from the moon surface. The transfer time from the resonant Cycler orbit to the earth is 0.5 day, and the height of the endpoint of the Lambert transfer orbit from the surface of the earth is 100 km. Assuming that the earth-to-resonant Cycler orbit transfer time is 0.5 days, the Lambert transfer orbit starting point is 100km from the earth's surface. The transfer time from the resonant Cycler orbit to the moon is 1 day, and the height from the end point of the Lambert transfer orbit to the surface of the moon is 100 km.

Using a Gaussian-global variable composite algorithm to pair SM_BCMProcessing to obtain a first set of Lambert transfer tracksThe corresponding Lambert track simulation results are shown in fig. 5. In the figure, under the earth-moon mass center inertial coordinate system O-XYZ, there are 4 broken lines, the first oneA Lambert transfer track from the Cycler track to the moon parking track; second dotted lineA Lambert transfer track from the moon parking track to the Cycler track; third dotted lineA Lambert transfer orbit from the earth mooring orbit to the Cycler orbit; fourth dotted lineFrom the Cycler orbit to the earth mooring orbitLambert transition tracks. Solid line is SM_BCMA track.

The classical two-body Lambert problem is a typical two-point boundary value problem, and a differential correction method is adopted to correct a first set of Lambert transfer tracksThe correction process is as follows:

according to L₁、L₂And θ can calculate the two-body Lambert transition tracks. The Lambert orbital transfer strategy gives a transfer orbit of position intersection, and the controlled quantity is an initial endpoint E₁Pulse velocity delta at position with directional angle noted β, differential correction algorithm will improve initial position velocity delta av₁And a transfer time t_fTo realize a terminal point E₂A meeting at a location. The convergence of the iterative process of the differential correction algorithm can be ensured by approaching the iterative initial value of the true value. For guiding a spacecraft P to a target position L₂(i.e., theoretical terminal position vector), each iteration integrates the orbit to L₂A, and t_fI.e. taken as the orbit integration time. Obviously, the initial endpoint E₁Track velocity v at a location₁₀Change of (Δ v)₁Will result in an orbital integration time t_fChange of (1), noted as Δ t_f. Examine the mth iteration (the previous iteration m times is recorded as m-1, the next iteration m times is recorded as m +1, then the iteration m times is recorded as the current time), and set the endpoint E₂The position vector at the position is first-order Taylor expanded to obtain:

$L(t_f+{\mathrm{\Δt}}_f,v_{10}+{\mathrm{\Δv}}_1)=L_2+\frac{\∂L_2}{\∂v_{10}}\×{\mathrm{\Δv}}_1+v_{20}\×{\mathrm{\Δt}}_f---\left(7\right)$

L(t_f+Δt_f,v₁₀+Δv₁) Representing the actual vector of the position of the spacecraft P terminal in an ellipse transferring orbit E relative to the transferring time₁A function between the initial velocities at the points;

representing terminal position vectors and elliptical transfer orbits E₁The partial derivative between the initial velocities at a point, reduces to

In the present invention, the speed correction amount Δ v₁Should be such that:

L(t_f+Δt_f,v₁₀+Δv₁)＝L₂(8)

in finishing formulae (7) and (8), one can obtain:

${\mathrm{\δL}}_2=L_2-\stackrel{\‾}{L_2}=-\frac{\∂L_2}{\∂v_{10}}\×{\mathrm{\Δv}}_1-v_{20}\×{\mathrm{\Δt}}_f---\left(9\right)$

representing the actual terminal position vector of the spacecraft P.

Is not to takeAnd L₂With the same abscissa component, the coordinate system O is rotated in the sun_S-X_SY_SZ_SNext, the terminal position vector and the ellipse transfer orbit E are recorded₁Partial derivative between initial velocities at pointsExpansion (9) yields:

$\left[\begin{array}{c}0\\ \mathrm{\δ}y\\ \mathrm{\δ}z\end{array}\right]=-\left[\begin{array}{cccc}{M}_{11}& M_{12}& M_{13}& v_{20}^x\\ M_{21}& M_{22}& M_{23}& v_{20}^y\\ M_{31}& M_{32}& M_{33}& v_{20}^z\end{array}\right]\×\left[\begin{array}{c}\mathrm{\Δ}{v}_1^x\\ {\mathrm{\Δv}}_1^y\\ {\mathrm{\Δv}}_1^z\\ {\mathrm{\Δt}}_f\end{array}\right]---\left(10\right)$

y denotes a daily rotation coordinate system O_S-X_SY_SZ_SLower Y_SA difference in position of the shaft;

z represents a rotational coordinate system of the sun and the earth_S-X_SY_SZ_SLower Z_SA difference in position of the shaft;

M₁₁representing terminal position vectors and elliptical transfer orbits E₁A first row, a first column element of a partial derivative matrix between initial velocities at points;

M₁₂representing terminal position vectors and elliptical transfer orbits E₁A first row and a second column of elements of a partial derivative matrix between initial velocities at the points;

M₁₃representing terminal position vectors and elliptical transfer orbits E₁The first row, third column elements of the partial derivative matrix between the initial velocities at the points;

M₂₁representing terminal position vectors and ellipsesTransfer track E₁A second row and a first column element of a partial derivative matrix between initial velocities at the points;

M₂₂representing terminal position vectors and elliptical transfer orbits E₁A second row and a second column of elements of the partial derivative matrix between the initial velocities at the points;

M₂₃representing terminal position vectors and elliptical transfer orbits E₁Second row, third column element of partial derivative matrix between initial velocities at points;

M₃₁representing terminal position vectors and elliptical transfer orbits E₁A third row, first column element of the partial derivative matrix between the initial velocities at the points;

M₃₂representing terminal position vectors and elliptical transfer orbits E₁A third row, second column element of the partial derivative matrix between the initial velocities at the points;

M₃₃representing terminal position vectors and elliptical transfer orbits E₁The third row and column three elements of the partial derivative matrix between the initial velocities at the points;

indicating the terminal position E₂End speed at X_SA velocity component of the shaft;

indicating the terminal position E₂At ending speed of Y_SA velocity component of the shaft;

indicating the terminal position E₂At ending speed of Z_SA velocity component of the shaft;

represents the initial endpoint E₁Initial velocity increment at position X_SA velocity component of the shaft;

represents the initial endpoint E₁Initial velocity increment at position Y_SA velocity component of the shaft;

represents the initial endpoint E₁Initial velocity increment at position Z_SA velocity component of the shaft;

Δt_frepresents the initial endpoint E₁Track velocity v at a location₁₀Delta of Δ v₁Resulting in a change in the track integration time.

Order toThen there are:

$\left[\begin{array}{c}{\mathrm{\Δv}}_1^x\\ {\mathrm{\Δv}}_1^y\\ {\mathrm{\Δv}}_1^z\\ {\mathrm{\Δt}}_f\end{array}\right]=-{\left(C^{T}C\right)}^{-1}{C}^T\left[\begin{array}{c}0\\ \mathrm{\δ}y\\ \mathrm{\δ}z\end{array}\right]---\left(11\right)$

C^Trepresenting the inverse of the matrix C.

Therefore, the speed increment of the current time m is characterized as follows relative to the speed of the previous time m-1:

${\mathrm{\Δv}}_1=\left[\begin{array}{ccc}{\mathrm{\Δv}}_1^x& {\mathrm{\Δv}}_1^y& {\mathrm{\Δv}}_1^z\end{array}\right]---\left(12\right)$

and repeating iteration after the accuracy is set until the accuracy is met, wherein the accuracy is set to be 1 × 10 for the transfer orbit of the moon berthing orbit taking off to the resonance type circulator orbit^－8The iteration results are shown in the table below.

Iteration 1	2.0620622902499468×100
		Iteration 2	2.7641929155045797×10-1
Iteration 3	7.6981072705594636×10-2
		Iteration 4	8.3862124467837824×10-3
Iteration 5	1.7139690673584045×10-4
		Iteration 6	1.0442627289698014×10-7

And the set precision is achieved through 6 iterations.

The accuracy of the transfer orbit from the resonance type circulator orbit to the moon parking orbit was set to 1 × 10^－8The iteration results are shown in the table below.

Iteration 1	2.5421929920011647×100
		Iteration 2	1.2546155043467997×10-1
Iteration 3	6.5637444323567223×10-2
		Iteration 4	5.3832424435867889×10-3
Iteration 5	2.2263784592810393×10-4
		Iteration 6	2.0143221442658014×10-7

And the set precision is achieved through 6 iterations.

The set precision is 1 × 10 aiming at the transfer orbit from the earth mooring orbit to the resonance type circulator orbit^－8The iteration results are shown in the table below.

Iteration 1	8.5792468327776916×10-1
		Iteration 2	7.9205098920078354×10-1
Iteration 3	1.4596634587845803×10-1
		Iteration 4	2.5378649526028702×10-3
Iteration 5	7.2114508959175548×10-7

And 5 times of iteration is performed to achieve the set precision.

The accuracy of the transfer orbit from the resonance type Cycler orbit to the earth mooring orbit is set to 1 × 10^－8The iteration results are shown in the table below.

And 7 times of iteration is performed to achieve the set precision.

In the present invention, a differential correction algorithm pair is adoptedCorrecting to obtain a first set of corrected Lambert transfer tracksAs shown in fig. 5A and 5B: in the figure, the correction is a thick solid line. The flying condition of the spacecraft on the actual transfer orbit is closer to the thick solid line and more accords with the motion rule. Compare FIG. 5 with FIG. 5A, whereinIs composed ofThe correction of (2) is carried out,is composed ofAnd (4) correcting. Compare FIG. 5 with FIG. 5B, whereinIs composed ofThe correction of (2) is carried out,is composed ofAnd (4) correcting.

Analysis of the intersection of the spacecraft take-off from the ground to the Cycler orbit

Emission window analysis was performed first:

the intersection angle of the white and red (the equatorial plane of the earth) changes in a weekly and monthly way, and the change range from 3 and 21 days in 2005 to 3 and 21 days in 2006 is [0 degrees and eta ], wherein the eta value changes from 28.4 degrees to 28.7 degrees. The launch latitude of the Wenchang is 27.9 degrees, the launch latitude of the Wenchang is 19 degrees, both the launch latitude and the launch latitude have the launch opportunity to directly launch the space station into the white road surface, and the window of the Wenchang for direct track entry is smaller than that of the Wenchang; as shown in fig. 6, the abscissa is time, and the ordinate is the intersection angle of white and red (the earth equatorial plane). For Wenchang, the emission window is about 7 days wide with an interval of about 7 days; for West Chang, the emission window is about 2 days wide and about 13 days apart. Within the emission window, the appropriate emission timing (1 time per day) should be chosen to capture the ascending crossover point into the white channel.

The launching latitude of West Chang refers to the launch site of West Chang who "Dragon Water" published in 2007-sketch of four satellite launching centers (1) in China, page 17, authored by West Chang.

The launching latitude of Wenchang refers to the thunder and lightning environment analysis of launching site of Hainan Wenchang rocket published in 2012, page 183, high YI by authors, etc.

Then a rendezvous window analysis is performed:

the space station was transferred from the parking track (height 100km) to the resonant Cycler track, and fuel consumption and transfer time were analyzed in a huffman transfer mode or a Lambert transfer mode, as shown in fig. 6A and 6B. The Hotman shift for the most fuel efficient in FIG. 6A requires 3192m/s for the velocity increment and 7 days for the shift time. The required velocity increase for the Hotman transfer of the top-cost fuel in FIG. 6B is 4093.5m/s, with a transfer time of 66.5 min.

Since fig. 6A and 6B arrive at apogee approximately equal in time (about 7 days), while fig. 6A is more fuel efficient, it will be taken as the nominal meeting trajectory: marking the nominal intersection as a launching track of a day 1, and if the days 2, 3, 4, 5, 6 and 7 are selected for launching because of the delay of the nominal intersection launching, sequentially increasing the required fuel from 3192km/s to 4093.5 m/s; marking the nominal intersection as a 7 th-day emission track, and optionally counting from several days to-6, -5, -4, -3, -2 and-1 days for emission, wherein the required fuel is increased to 4093.5m/s from 3192km/s in sequence; and marking the nominal intersection as the emission track of 1-7 days, and increasing the required fuel from 3192km/s to 4093.5m/s in sequence. The emission locus for 7 consecutive days is shown in fig. 6C and 6D, the emission locus under the Earth-moon mass center inertial coordinate system O-XYZ is shown in fig. 6C, and the emission locus under the Earth-moon mass center rotational coordinate system O-XYZ is shown in fig. 6D.

(II) analysis of the intersection of the spacecraft from the lunar takeoff to the resonant Cycler orbit

Emission window analysis was performed first:

the intersection angle of the white and red (the equatorial plane of the moon) varies in the weekend, and the range of variation from 3/21/2005 to 3/21/2006 is [0 ° ], in which the value varies from 6.64 ° to 6.86 °, and as shown in fig. 7, the abscissa is time and the ordinate is the intersection angle of the white and red (the equatorial plane of the moon). The variation period was about 13.6 days. For a landing point with the latitude higher than 6.64 degrees, the plane of the orbit needs to be adjusted after the lunar takeoff, and the landing point can enter the orbit of the resonant circulator. Taking a south-north polar landing point as an example, the lunar surface takes off and enters a 100km circular lunar polar orbit; the resonant Cycler orbital capture process is as follows: lifting the moon point to 20000km, wherein the required speed increment is about 577 m/s; the moon point was tilt corrected and the required velocity increment was approximately 287.3 m/s. The lifting of the moon point saves part of the energy for the next resonant Cycler orbital transaction. The revolution period and the rotation period of the moon are the same, and only two times of each month can capture the ascension point into the white road surface.

Then a rendezvous window analysis is performed:

and searching a meeting window with the speed increment less than or equal to 1500m/s by taking the meeting time and the phase of the resonant Cycler orbit as variables and taking the height of the lunar parking orbit of 100km as a constraint condition. The geographical longitude depends on the ascent crossing right ascension, and the search result according to the "height of parking track 100 km" as the constraint condition is shown in fig. 7A, where in fig. 7A, points T1, T2, T3 and T4 represent the optimal meeting conditions that meet the search condition. The geographical longitudes as shown in fig. 7B may be distributed in [ -170 °, -135 ° ] [ -43 °, -5 ° ] [ -45 °,108 ° ] [150 °,165 ° ]. The landing points in the interval are expected to directly enter the white road surface after the landing points take off through the lunar surface; and the landing points outside the interval need to wait for half a month for the ascending crossing right ascension capture. Therefore, if a landing point in the designated area is selected, the moon surface enters the white road surface after taking off: the phase of resonant circulator orbit is [0.05,0.25 ]. sup.U [0.8,1], so that it can take off and meet at any time.

In the invention, a Gaussian-global variable composite algorithm is adopted to carry out SM_LProcessing to obtain a second set of Lambert transfer tracksFor SM_LIs processed by the method and the pair SM_BCMThe treatment method of (2) is the same. For the purpose of distinguishing between the descriptions, the letter H is added to the cited relational expression for distinction.

According to Lambert's law of flight time, the transfer time between any two points on the elliptical transfer orbit and the sum (HL) of the major semi-axis Hra and the radius of the two points of the elliptical transfer orbit₁+HL₂) And center angle H θ, then expressed as:

Ht_f＝W(Hra,(HL₁+HL₂),HR_P) (13)

if HL₁And HL₂The sum is constant, the major semi-axis Hra is constant, and the initial end point E₁And terminal E₂Distance between (HR)_PIs constant, then from the initial end point E₁To terminal E₂Time of flight Ht_fIs also a constant. The determination of the elliptical transfer orbit and the choice of the speed between the two points are the key to the Lambert problem, which is described as a gaussian problem as follows: tracking spacecraft P departure position vector (HL)₁) And velocity vector (Hv)₁) Spacecraft P terminal position vector (HL)₂) And velocity vector (Hv)₂) Flight transition time of Ht_fWith spacecraft P in ellipseTransfer track E₁The initial velocity at a point is Hv₁₀The spacecraft P in the elliptical transfer orbit E₂The ending velocity at the point is Hv₂₀. The goal of the Gaussian problem is to solve for the initial position velocity increment Δ Hv₁Velocity increment delta Hv from terminal position₂And the magnitude of the applied pulse thrust is determined.

In the present invention, the gaussian problem can be solved by the transcendental system of equations as follows:

HL₂＝k×HL₁+g×Hv₁₀(14)

${\mathrm{Hv}}_{20}=\stackrel{\·}{k}\×{\mathrm{HL}}_1+\stackrel{\·}{g}\×{\mathrm{Hv}}_{10}---\left(15\right)$

lagrange coefficient of one

Lagrange coefficient of two

Lagrange coefficient of three

Lagrange coefficient of fourMu is the mass ratio of the secondary celestial body to the primary celestial body, and in a double circle model BCM, the value is 0.000003003; h is a variable on the elliptical transfer orbit of the spacecraft P at different positions.

In the formula (14), the formula (15), HL is known₁And Hv₁₀Or HL, or₂And Hv₂₀The elliptical transfer orbit followed by the spacecraft P can be determined. It is clear that once the Lagrangian coefficients k, g are determined,The gaussian problem can be readily solved.

In the invention, the formula (14) and the formula (15) are solved by adopting a global variable algorithm.

Assuming that the transfer time from the earth to the co-homed orbit is 1 day, the height of the starting point of the transfer orbit from the earth surface is 100 km. The transfer time from the co-homed orbit to the moon is 4 days, and the height of the transfer orbit terminal point to the moon surface is 100 km. The corresponding homographic Lambert track simulation results are shown in fig. 7C according to a two-body Lambert problem solution.

Using differential correction algorithm pairCorrecting to obtain a second set of corrected Lambert transfer tracksThe second set of corrected transfer orbit simulation results are shown in fig. 7D.

According to HL₁、HL₂And H θ can calculate the two-body Lambert transition track. The Lambert orbital transfer strategy gives a transfer orbit of position intersection, and the controlled quantity is an initial endpoint E₁The pulse velocity increment at position, its direction angle is denoted as H β, the differential correction algorithm will improve the initial position velocity increment Δ Hv₁And transfer time Ht_fTo realize a terminal point E₂A meeting at a location. The iterative initial value close to the true value can ensure the differential correction algorithmAnd (5) convergence of the iteration process. For guiding a spacecraft P to a target position L₂(i.e., theoretical terminal position vector), each iteration integrates the orbit into HL₂And Ht_fI.e. taken as the orbit integration time. Obviously, the initial endpoint E₁Track velocity Hv at location₁₀Change of (a) Δ Hv₁Will result in an orbital integration time Ht_fChange of (d), noted as Δ Ht_f. Examine the mth iteration (the previous iteration m times is recorded as m-1, the next iteration m times is recorded as m +1, then the iteration m times is recorded as the current time), and set the endpoint E₂The position vector at the position is first-order Taylor expanded to obtain:

$HL({\mathrm{Ht}}_f+{\mathrm{\ΔHt}}_f,{\mathrm{Hv}}_{10}+{\mathrm{\ΔHv}}_1)={\mathrm{HL}}_2+\frac{\∂{\mathrm{HL}}_2}{\∂{\mathrm{Hv}}_{10}}\×{\mathrm{\ΔHv}}_1+{\mathrm{Hv}}_{20}\×{\mathrm{\ΔHt}}_f---\left(16\right)$

HL(Ht_f+ΔHt_f,Hv₁₀+ΔHv₁) Representing space flightThe actual vector of the position of the terminal of the device P is transferred in an ellipse on the orbit E with respect to the transfer time₁A function between the initial velocities at the points;

representing terminal position vectors and elliptical transfer orbits E₁The partial derivative between the initial velocities at a point, reduces to

In the present invention, the speed correction amount Δ Hv₁Should be such that:

HL(Ht_f+ΔHt_f,Hv₁₀+ΔHv₁)＝HL₂(17)

in accordance with formulae (16) and (17), the following can be obtained:

${\mathrm{\δHL}}_2={\mathrm{HL}}_2-\stackrel{\‾}{{\mathrm{HL}}_2}=\frac{\∂{\mathrm{HL}}_2}{\∂{\mathrm{HV}}_{10}}\×{\mathrm{\ΔHv}}_1-{\mathrm{Hv}}_{20}\×{\mathrm{\ΔHt}}_f---\left(18\right)$

representing the actual terminal position vector of the spacecraft P.

Is not to takeAnd L₂With the same abscissa component, the coordinate system O is rotated in the sun_S-X_SY_SZ_SNext, the terminal position vector and the ellipse transfer orbit E are recorded₁Partial derivative between initial velocities at pointsExpansion (18) yields:

$\left[\begin{array}{c}0\\ \mathrm{\δ}Hy\\ \mathrm{\δ}Hz\end{array}\right]=-\left[\begin{array}{cccc}{\mathrm{HM}}_{11}& {\mathrm{HM}}_{12}& {\mathrm{HM}}_{13}& {\mathrm{Hv}}_{20}^x\\ {\mathrm{HM}}_{21}& {\mathrm{HM}}_{22}& {\mathrm{HM}}_{23}& {\mathrm{Hv}}_{20}^y\\ {\mathrm{HM}}_{31}& {\mathrm{HM}}_{32}& {\mathrm{HM}}_{33}& {\mathrm{Hv}}_{20}^z\end{array}\right]\×\left[\begin{array}{c}{\mathrm{\ΔHv}}_1^x\\ {\mathrm{\ΔHv}}_1^y\\ {\mathrm{\ΔHv}}_1^z\\ {\mathrm{\ΔHt}}_f\end{array}\right]---\left(19\right)$

hy represents a daily-ground rotation coordinate system O_S-X_SY_SZ_SLower Y_SA difference in position of the shaft;

hz denotes the sun-ground rotation coordinate system O_S-X_SY_SZ_SLower Z_SA difference in position of the shaft;

HM₁₁representing terminal position vectors and elliptical transfer orbits E₁A first row, a first column element of a partial derivative matrix between initial velocities at points;

HM₁₂representing terminal position vectors and elliptical transfer orbits E₁A first row and a second column of elements of a partial derivative matrix between initial velocities at the points;

HM₁₃representing terminal position vectors and elliptical transfer orbits E₁The first row, third column elements of the partial derivative matrix between the initial velocities at the points;

HM₂₁representing terminal position vectors and elliptical transfer orbits E₁Of the initial velocity at the pointSecond row first column elements of the partial derivative matrix in between;

HM₂₂representing terminal position vectors and elliptical transfer orbits E₁A second row and a second column of elements of the partial derivative matrix between the initial velocities at the points;

HM₂₃representing terminal position vectors and elliptical transfer orbits E₁Second row, third column element of partial derivative matrix between initial velocities at points;

HM₃₁representing terminal position vectors and elliptical transfer orbits E₁A third row, first column element of the partial derivative matrix between the initial velocities at the points;

HM₃₂representing terminal position vectors and elliptical transfer orbits E₁A third row, second column element of the partial derivative matrix between the initial velocities at the points;

HM₃₃representing terminal position vectors and elliptical transfer orbits E₁The third row and column three elements of the partial derivative matrix between the initial velocities at the points;

indicating the terminal position E₂End speed at X_SA velocity component of the shaft;

indicating the terminal position E₂At ending speed of Y_SA velocity component of the shaft;

indicating the terminal position E₂At ending speed of Z_SA velocity component of the shaft;

represents the initial endpoint E₁Initial velocity increment at position X_SA velocity component of the shaft;

represents the initial endpoint E₁Initial velocity increment at position Y_SA velocity component of the shaft;

represents the initial endpoint E₁Initial velocity increment at position Z_SA velocity component of the shaft;

ΔHt_frepresents the initial endpoint E₁Track velocity v at a location₁₀Delta of Δ v₁Resulting in a change in the track integration time.

Order toThen there are:

$\left[\begin{array}{c}{\mathrm{\ΔHv}}_1^x\\ {\mathrm{\ΔHv}}_1^y\\ {\mathrm{\ΔHv}}_1^z\\ {\mathrm{\ΔHt}}_f\end{array}\right]=-{\left({\mathrm{HC}}^{T}HC\right)}^{-1}{\mathrm{HC}}^T\left[\begin{array}{c}0\\ \mathrm{\δ}Hy\\ \mathrm{\δ}Hz\end{array}\right]---\left(20\right)$

HC^Trepresenting the inversion of the matrix HC.

Therefore, the speed increment of the current time m is characterized as follows relative to the speed of the previous time m-1:

${\mathrm{\ΔHv}}_1=\left[\begin{array}{ccc}{\mathrm{\ΔHv}}_1^x& {\mathrm{\ΔHv}}_1^y& {\mathrm{\ΔHv}}_1^z\end{array}\right]---\left(21\right)$

and repeating iteration after the accuracy is set until the accuracy is met, wherein the accuracy is set to be 1 × 10 for the transfer orbit of the moon berthing orbit taking off to the resonance type circulator orbit^－15The iteration results are shown in the table below.

Iteration 1	2.8639968709823611×100
		Iteration 2	1.2827923060663518×100
Iteration 3	5.6581823976188161×10-1
		Iteration 4	2.1973874406340671×10-1
Iteration 5	1.5468617860521472×10-1
		Iteration 6	4.9764011654783211×10-2
Iteration 7	5.2642496344864611×10-3
		Iteration 8	2.0846898885215647×10-5
Iteration 9	2.2579833858870749×10-10
		Iteration 10	6.5486562934038266×10-14

And after 10 iterations, the set precision is achieved.

The precision is set to be 1 × 10 aiming at the transfer orbit from the earth mooring orbit to the co-homed Cycler orbit^-11The iteration results are shown in the table below.

Iteration 1	4.0437979717451373×100
		Iteration 2	2.1369334069815822×100
Iteration 3	9.8469011235887069×10-2
		Iteration 4	5.1247925851644960×10-4
Iteration 5	5.5164314416937229×10-8
		Iteration 6	3.3282881768155430×10-12

And the set precision is achieved through 6 iterations.

In the present invention, the comparison is based on a resonant Cycler orbit SM_BCMAnd co-lodging type Cycler orbit SM_LTime of flight (including orbital transfer time), fuel consumption of the earth-moon transfer mission. Fuel consumption was characterized by the total velocity momentum of two Lambert rail changes per task.

1) Resonant Cycler orbital protocol:

the total flight time from the earth mooring trajectory to the moon mooring trajectory was 6.014 days.

The total velocity impulse of the two Lambert orbits during this period was 4.626 km/s.

2) The same-host Cycler orbit scheme:

the total flight time from the earth mooring trajectory to the moon mooring trajectory was 20.695 days.

The total velocity impulse of the two Lambert orbits during this period was 4.450 km/s.

By comparing the total flight time and the total fuel consumption, it can be seen that the earth-moon transfer scheme based on the Cycler orbit has longer total flight time and less total fuel consumption. This means that based on the Cycler orbitThe concept of the march transfer scenario is feasible. Both types of Cycler orbits have good periodicity. The period of the resonant circulator orbit can be set, and the spacecraft can regularly transmit materials and personnel to and from the earth and the moon in a planned way according to the task requirement. The period of the cycle of the homeopathic Cycler orbit can also be determined and passed through the Earth-moon line LL₁When the round-trip transport task of the earth and the moon is finished, the earth and the moon can be tied LL₁And replacing and overhauling the scientific equipment arranged at the point. The earth-moon transfer scheme based on the two types of Cycler orbits has longer flight time and less total fuel consumption. The advantage of good periodicity is reasonably utilized, a moon-earth round-trip public traffic system can be designed, and the system has great significance for deep space exploration of a low-thrust spacecraft and material transmission of the moon-earth system.

Claims (7)

1. A earth-moon round trip task simulation system based on a Cycler orbit is characterized in that: the simulation system comprises a resonant Cycler orbit model building module (10), a co-homed Cycler orbit model building module (20), a multi-target method orbit correction module (30) and a Lambert transfer orbit acquisition module (50);

the first aspect of the method for constructing the resonant Cycler orbit model (10) is that M is constructed according to a circular restrictive three-body problem model CR3BP_CR3BPA kinetic model; in the second aspect, the double-circle model BCM is adopted to match the M_CR3BPOptimizing the dynamic model to obtain M_BCMA kinetic model;

the construction of the co-homed Cycler orbit model (20) is based on the double circle model BCM to construct M_LA kinetic model;

the multiple targeting method track correction module (30) adopts the multiple targeting method to respectively correct M_BCMCorrecting to obtain a corrected resonance Cycler orbit dynamics model DM_BCM(ii) a Second aspect to DM_BCMMethod for obtaining resonant Cycler orbit SM by adopting fourth-order Runge Kutta method_BCM(ii) a In the third aspect, multiple targeting methods are used to target M separately_LCorrecting to obtain a corrected compatible Cycler orbit dynamics model DM_L(ii) a Fourth aspect to DM_LObtaining co-sited Cycler orbit SM by adopting four-step Runge Kutta method_L；

The first aspect of the Lambert transition track acquisition module (50) employs a Gaussian-global variable composite algorithm to pair SM_BCMProcessing to obtain a first set of Lambert transfer tracksThe second aspect adopts a differential correction algorithm pairCorrecting to obtain a first set of corrected Lambert transfer tracksThird aspect is to measure SM according to variation rule of emission point latitude and white-red intersection angle_BCMAnalyzing a launching window and an intersection window to obtain the launching and orbit entering time of the spacecraft; the fourth aspect adopts a Gaussian-global variable composite algorithm to SM_LProcessing to obtain a second set of Lambert transfer tracksThe fifth aspect adopts a differential correction algorithm pairTo carry out repairPositive to obtain a second set of corrected Lambert transition tracksThe sixth aspect can learn the SM by comparing the fuel consumption amount and the total flight time of the earth-moon round trip mission_BCMAnd SM_LThe respective advantages of the orbits provide optimized design indexes for the earth-moon round trip task of the deep space exploration of the low-thrust spacecraft.

2. The Cycler orbital based ground-to-month round trip task simulation system of claim 1, wherein: in the Lambert transfer track acquisition module (50), according to L₁、L₂And θ can calculate the two-body Lambert transition track; the Lambert orbital transfer strategy gives a transfer orbit of position intersection, and the controlled quantity is an initial endpoint E₁Pulse velocity delta at position with directional angle noted β, differential correction algorithm will improve initial position velocity delta av₁And a transfer time t_fTo realize a terminal point E₂A meeting at a location; for guiding a spacecraft P to a target position L₂Each iteration integrates the orbit to L₂A, and t_fTaking the orbit integration time; obviously, the initial endpoint E₁Track velocity v at a location₁₀Change of (Δ v)₁Will result in an orbital integration time t_fChange of (1), noted as Δ t_f(ii) a Examine the mth iteration, for end point E₂The position vector at the position is first-order Taylor expanded to obtainRepresenting the actual vector of the position of the spacecraft P terminal in an ellipse transferring orbit E relative to the transferring time₁A function between the initial velocities at the points;representing terminal position vectors and elliptical transfer orbits E₁The partial derivative between the initial velocities at a point, reduces toL₁As an initial end point E of the spacecraft P₁A position vector of (a); l is₂Terminal E for a spacecraft P₂A position vector of (a); theta is a transfer angle;

speed correction amount Deltav₁Should be such that L (t)_f+△t_f,v₁₀+△v₁)＝L₂Is established, then

Represents the actual terminal position vector of the spacecraft P;

getAnd L₂With the same abscissa component, the coordinate system O is rotated in the sun_S-X_SY_SZ_SNext, the terminal position vector and the ellipse transfer orbit E are recorded₁Partial derivative between initial velocities at pointsThen there is

y denotes a daily rotation coordinate system O_S-X_SY_SZ_SLower Y_SA difference in position of the shaft;

z represents a rotational coordinate system of the sun and the earth_S-X_SY_SZ_SLower Z_SA difference in position of the shaft;

M₁₁representing terminal position vectors and elliptical transfer orbits E₁A first row, a first column element of a partial derivative matrix between initial velocities at points;

M₁₂representing terminal position vectors and elliptical transfer orbits E₁A first row and a second column of elements of a partial derivative matrix between initial velocities at the points;

M₁₃representing terminal position vectors and elliptical transfer orbits E₁The first row, third column elements of the partial derivative matrix between the initial velocities at the points;

M₂₁representing terminal position vectors and elliptical transfer orbits E₁A second row and a first column element of a partial derivative matrix between initial velocities at the points;

M₂₂representing terminal position vectors and elliptical transfer orbits E₁A second row and a second column of elements of the partial derivative matrix between the initial velocities at the points;

M₂₃representing terminal position vectors and elliptical transfer orbits E₁Second row, third column element of partial derivative matrix between initial velocities at points;

M₃₁representing terminal position vectors and elliptical transfer orbits E₁A third row, first column element of the partial derivative matrix between the initial velocities at the points;

M₃₂representing terminal position vectors and elliptical transfer orbits E₁A third row, second column element of the partial derivative matrix between the initial velocities at the points;

M₃₃representing terminal position vectors and elliptical transfer orbits E₁The third row and column three elements of the partial derivative matrix between the initial velocities at the points;

indicating the terminal position E₂End speed at X_SA velocity component of the shaft;

indicating the terminal position E₂At ending speed of Y_SA velocity component of the shaft;

indicating the terminal position E₂At ending speed of Z_SA velocity component of the shaft;

represents the initial endpoint E₁Initial velocity increment at position X_SA velocity component of the shaft;

represents the initial endpoint E₁Initial velocity increment at position Y_SA velocity component of the shaft;

represents the initial endpoint E₁Initial velocity increment at position Z_SA velocity component of the shaft;

△t_frepresents the initial endpoint E₁Track velocity v at a location₁₀Increment △ v₁The resulting change in the track integration time;

order toThen there isC^TRepresents the inversion of the matrix C; therefore, the velocity increment of the current time m is characterized as the velocity increment of the current time m relative to the velocity of the previous time m-1

3. The Cycler orbital based ground-to-month round trip task simulation system of claim 1, wherein: obtaining a resonant Cycler orbit dynamics model according to a circular restrictive three-body problem model CR3BPWherein,

represents the partial derivative in the x-direction;

u represents the potential function of the spacecraft P, andμ₁represents the mass ratio of the moon to the earth, generally taking the value of 0.01215; r₁Representing space P to earth P₁Distance of (A), R₂Representing spacecraft P to moon P₂The distance of (d);

represents the partial derivative in the y-direction;

represents the partial derivative in the z direction;

VA^xrepresents the velocity of the spacecraft P in the x direction;

VA^yrepresents the velocity of the spacecraft P in the y direction;

HA^xrepresents the acceleration of the spacecraft P in the x direction;

HA^yrepresents the acceleration of the spacecraft P in the y direction;

HA^zrepresenting the acceleration of the spacecraft P in the z direction.

4. The Cycler orbital based ground-to-month round trip task simulation system of claim 1, wherein: obtaining a co-homed Cycler orbit dynamics model according to a double-circle model BCMWherein,

in the double circle model BCM, a day-ground rotation coordinate system O is used_S-X_SY_SZ_SThe sun is P₃The spacecraft P rotates the coordinate system O in the sun and the earth_S-X_SY_SZ_SThe lower position is denoted as (xs)_P,ys_P,zs_P) (ii) a Moon P₂In the sun-ground rotation coordinate system O_S-X_SY_SZ_SThe lower position is notedr_PSRepresents the dimensionless distance of the spacecraft P from the sun, anr_PERepresents the dimensionless distance of the spacecraft P from the earth, anr_PMRepresents the dimensionless distance of the spacecraft P from the moon, anμ₂Representing the mass ratio of the earth to the sun; m is_MRepresenting the nondimensionalized quality of the moon; by the centroid of the sun and earth_SIs an origin, and the sun-ground connecting line is X_SAxis, pointing towards the earth, X_SThe positive direction of the axis; y is_SThe axis being perpendicular to the plane of motion of the sun and earth, and Z being determined according to the right-hand rule_SThe direction of the axis;

indicates M is_LSpacecraft under model P at X_SA speed in a direction;

indicates M is_LSpacecraft under model P inY_SA speed in a direction;

indicates M is_LSpacecraft under model P at X_SAcceleration in a direction;

indicates M is_LSpacecraft under model P in Y_SAcceleration in a direction;

indicates M is_LSpacecraft under model P in Z_SAcceleration in the direction.

5. The Cycler orbital based ground-to-month round trip task simulation system of claim 1, wherein: the method comprises the steps that the latitude of a launching point and the intersection angle change rule of white and red are utilized, a space station is transferred to a resonant Cycler track from a parking track, a Hoeman transfer mode or a Lambert transfer mode is adopted to analyze fuel consumption and transfer time, the needed speed increment of Hoeman transfer of the most-economical fuel is 3192m/s, and the transfer time is 7 days; the increment of the speed required by the Hulman transfer of the most-cost fuel is 4093.5m/s, and the transfer time is 66.5 min; taking Hotman transfer of the most fuel-saving fuel as a nominal intersection track, and if the 2 nd, 3 rd, 4 th, 5 th, 6 th and 7 th days are selected for emission because of the delay of the nominal intersection emission, sequentially increasing the required fuel from 3192km/s to 4093.5 m/s; and marking a nominal intersection as a 7 th-day emission track, or selecting to carry out emission from counting the number of the reciprocal days, namely-6, -5, -4, -3, -2 and-1 days, and increasing the required fuel from 3192km/s to 4093.5m/s in sequence.

6. The Cycler orbital based ground-to-month round trip task simulation system of claim 1, wherein: obtaining by using the latitude of a transmitting point and a change rule of a white-red intersection angle, taking intersection time and a resonant Cycler orbit phase as variables, taking the height of a lunar parking orbit of 100km as a constraint condition, and searching an intersection window with the speed increment less than or equal to 1500 m/s; the geographic longitude depends on the ascent crossing right ascension, and the optimal meeting condition can be searched according to the height of the parking track of 100km as a constraint condition.

7. The Cycler orbital based ground-to-month round trip task simulation system of claim 1, wherein: the comparison is based on a resonant Cycler track SM_BCMAnd co-lodging type Cycler orbit SM_LThe flight time and fuel consumption of the earth-moon transfer task; representing the fuel consumption by the total speed impulse of two Lambert rail changes of each task;

the resonance type Cycler orbit: the total flight time from the earth mooring trajectory to the moon mooring trajectory was 6.014 days; the total velocity impulse of two Lambert orbital changes in the period is 4.626 km/s;

the homeopathic Cycler orbit: the total flight time from the earth mooring trajectory to the moon mooring trajectory was 20.695 days; the total velocity impulse of the two Lambert orbits during this period was 4.450 km/s.