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

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

Technical field

The present invention relates to the task simulation come and gone between a kind of earth and the moon, more particularly, refer to a kind of based on The task simulation system that the ground moon of Cycler track comes and goes.

Background technology

First Chinese moon exploration project is by moon probing satellite, carrier rocket, launching site, observing and controlling and Ground Application etc. five Big system composition, China's lunar exploration satellite engineering also has five large-engineering targets: one is to develop and launch first lunar exploration satellite of China； Two is tentatively to grasp lunar orbiting exploration basic fundamental；Three is to carry out lunar science detection first；Four is Primary Construction moon exploration boat It engineering system；Five is to accumulate experience for moon exploration successive projects.Key technology for this Gonna breakthrough moon probing satellite；Just Step sets up the big system of deep space exploration program of China；Checking every key technology such as payload and data interpretation；Preliminary foundation China's survey of deep space technology develops system；Cultivate the corresponding talent team.

First Chinese moon exploration project four big science task is:

One, moonscape 3 D stereoscopic image is obtained, the fine essential structure dividing moonscape and geomorphic unit, carry out The research of moonscape impact crater form, size, distribution, density etc., for division and the evolution in early days of terrestrial planet surface age Historic survey provides master data, and preferably provides basic data etc. for lunar surface soft landing district's addressing and position, lunar base.

Two, analyzing moonscape useful element content and the characteristic distributions of material type, mainly prospecting moonscape has The content of 14 kinds of elements such as the titanium of value of exploiting and utilizing, iron and distribution, draw the ball distribution map whole month of each element, lunar rock, Mineral and geology thematic map etc., find each element enrichment region at menology, the exploitation prospect of assessment moon mineral resources Deng.

Three, detection lunar soil thickness, i.e. utilizes microwave irradiation technology, obtains the thickness data of moonscape lunar soil, thus obtains To moonscape age and distribution thereof, and on this basis, content, resource distribution and the money of estimation nuclear fusion fuel used to generate electricity helium 3 Source amount etc..

Four, the detection earth is to the space environment of the moon.The moon and earth average distance are 380,000 kilometers, are in magnetic field of the earth The remote magnetic tail region in space, satellite is at this region detectable solar cosmic ray high energy particle and solar wind plasma body, and research is too The interaction of sun wind and the moon and magnetic field of the earth magnetic tail and the moon.

Cycler track refers to periodically travel to and fro between the earth and the moon, is diversion and the rail that do not stops near planet Road.Run on the aircraft on cyclic track can (several years even more than ten years) keep flying without track between planet for a long time Motor-driven (or having only to the least orbit maneuver), thus be considered as the one economical long-range mission mode saving energy. Cycler track scheme also can be sub-divided into whole circle track, half-turn track, general track etc. according to around day number of turns difference.Cycler Track can be divided into resonance type Cycler track and chummage type Cycler track.

Circular re stricted three body problem (Circular Restricted Three-Body Problem, CR3BP) describes Quality the most infinitesimal 3rd body fortune under the graviational interaction of two primary bodys moved in a circle around its public barycenter Dynamic.With reference in May, 2010, " Postgraduate School, Chinese Academy of Sciences Ph.D. Dissertation " of Li Mingtao, the phase of page 19 to the 21st page Close content.

Double round models (Bi-Circalar Model, BCM) are that one infinitely small mass body of research is in the moon (system for winding matter The heart operates on circular orbit), the gravitation effect of the system of the sun and the earth (operating on circular orbit around common barycenter) Under the basic model of the characteristics of motion.

Summary of the invention

The purpose of the present invention is to propose to the task simulation system come and gone by a kind of ground based on Cycler track moon, this system exists Resonance type Cycler track and chummage type Cycler track is produced under circular Restricted three-body model.Bias due to lunar orbit Rate and the impact of other celestial bodies (such as the sun, Jupiter etc.) gravitation so that the track set up under circular Restricted three-body model is with true Situation is different, and even due to the impact of perturbation, these tracks cannot keep constant.Therefore, with at circular Restricted three-body model The track of lower foundation is initial value, is optimized under double round models, utilize multiple shooting method to revise resonance type periodic orbit and Chummage type periodic orbit.Revised resonance type periodic orbit and chummage type periodic orbit can periodically meet with the earth, the moon, Create the launch window as target with the earth or the moon.For completing complete Earth-moon transfer orbit, analogue system of the present invention is also According to the classical solution of disome Lambert problem, generate the earth take off, " earth-cycle transfer orbit " of moon landing； " the cycle transfer orbit-moon " that the moon takes off, the earth lands；And be corrected under four body Models.Emulate in the present invention System calculates and compares Fuel Consumption and the flight time of both transfer orbits.

The present invention is the task simulation system come and gone by a kind of ground based on Cycler track moon, and this analogue system includes structure Build resonance Cycler model trajectory 10, build chummage Cycler model trajectory 20, multiple shooting method orbital exponent module 30 and Lambert transfer orbit acquisition module 50.

Described resonance Cycler model trajectory 10 first aspect that builds is according to circular re stricted three body problem MODEL C R3BP structure Build M_CR3BPKinetic model；Second aspect uses double circle Model B CM to described M_CR3BPKinetic model is optimized, and obtains M_BCMKinetic model.

Described structure chummage Cycler model trajectory 20 builds M according to double round Model B CM_LKinetic model.

Described multiple shooting method orbital exponent module 30 first aspect uses multiple shooting method respectively to M_BCMIt is modified, Obtain revised resonance Cycler dynamics of orbits model DM_BCM；Second aspect is to DM_BCMEmploying Fourth order Runge-Kutta obtains Resonance Cycler track SM_BCM；The third aspect uses multiple shooting method respectively to M_LIt is modified, obtains revised chummage Cycler dynamics of orbits model DM_L；Fourth aspect is to DM_LFourth order Runge-Kutta is used to obtain chummage Cycler track SM_L。

Described Lambert transfer orbit acquisition module 50 first aspect uses Gauss-global variable composite algorism to SM_BCM Process, obtain first set Lambert transfer orbitSecond aspect uses Differential correction algorithm pairIt is modified, obtains the revised Lambert of first set Transfer orbitThe third aspect becomes according to launch latitude and the white red angle of cut Law is to SM_BCMCarry out launch window and the analysis of intersection window, obtain Spacecraft Launch and the opportunity entered the orbit；Fourth aspect is adopted By Gauss-global variable composite algorism to SM_LProcess, obtain the second set Lambert transfer orbit 5th aspect uses differential correction algorithm pairIt is modified, obtains the second set revised Lambert transfer Track6th aspect by the Fuel Consumption of comparatively moon two-way mission, total flight time can Know SM_BCMAnd SM_LThe respective advantage of track, the ground moon two-way mission for the survey of deep space of low-thrust spacecraft provides optimization to set Meter index.

The advantage of analogue system of the present invention is:

1. the present invention produces resonance type Cycler track and chummage type Cycler track under circular Restricted three-body model. Be in order to revise existing resonance type Cycler track and chummage type Cycle track do not account for solar gravitation perturbation track is set The impact of meter.

2. present system considers solar gravitation perturbation when obtaining ground moon Cycler track, and uses multiple target practice Method obtains approximate period Cycler track, and precision increases compared with the cycle Cycler track that existing method calculates.

3. present system is taken off to the intersection window of Cycler track and spacecraft by the moon by ground analyzing spacecraft Face take off the intersection window to resonance Cycler track time, it is provided that a kind of iterative algorithm, solar gravitation perturbation feelings can considered Obtaining transition window under condition, precision increases compared with existing Lambert method result of calculation.

Accompanying drawing explanation

Fig. 1 is the structured flowchart of the task simulation system that present invention ground based on the Cycler track moon comes and goes.

Fig. 2 is ground moon barycenter inertial coodinate system O-XYZ and the coordinate system schematic diagram of ground moon barycenter rotating coordinate system O-xyz.

Fig. 2 A is day ground rotating coordinate system O_S-X_SY_SZ_SCoordinate system schematic diagram.

Fig. 3 is the multiple shooting method revised resonance Cycler dynamics of orbits model DM through the present invention_BCMFigure.

Fig. 3 A is the multiple shooting method revised chummage Cycler dynamics of orbits model DM through the present invention_LFigure.

Fig. 4 is the schematic diagram that under Lambert problem, 2 boundary values in astrodynamics obtain.

Fig. 5 is to process SM through Gauss of the present invention-global variable composite algorism_BCMAfter simulation result figure.

Fig. 5 A is the simulation result figure through differential correction algorithm revised first set correction Lambert transfer orbit.

Fig. 5 B is the simulation result figure revising Lambert transfer orbit through revised second set of differential correction algorithm.

Fig. 6 is the window schematic diagram directly entered the orbit in Xichang.

Fig. 6 A is the track schematic diagram of the required speed increment of the transfer of fuel saving Huo Man.

Fig. 6 B is the track schematic diagram of the required speed increment of transfer of expense fuel Huo Man.

Fig. 6 C is launching track schematic diagram under ground moon barycenter inertial coodinate system O-XYZ.

Fig. 6 D is launching track schematic diagram under ground moon barycenter rotating coordinate system O-xyz.

Fig. 7 is the window schematic diagram directly entered the orbit in Xichang that lunar surface takes off.

Fig. 7 A be parking orbit height be the Search Results schematic diagram under the constraints of 100km.

Fig. 7 B be geographic logitude be [-170 ° ,-135 °] ∪ [-43 ° ,-5 °] ∪ [45 °, 108 °] ∪ [150 °, 165 °] Enter the orbit window schematic diagram.

Fig. 7 C is the chummage Lambert track simulation result figure under disome Lambert problem.

Fig. 7 D is that Lambert transfer orbit revised by the second set after differential correction algorithmImitative True result.

Detailed description of the invention

Below in conjunction with accompanying drawing, the present invention is described in further detail.

Shown in Figure 1, the present invention is the task simulation system come and gone by a kind of ground based on Cycler track moon, this emulation System includes structure resonance Cycler model trajectory 10, builds chummage Cycler model trajectory 20, multiple shooting method track and repair Positive module 30 and Lambert transfer orbit acquisition module 50.

Described resonance Cycler model trajectory 10 first aspect that builds is according to circular re stricted three body problem MODEL C R3BP structure Build M_CR3BPKinetic model；Second aspect uses double circle Model B CM to described M_CR3BPKinetic model is optimized, and obtains M_BCMKinetic model.

Described structure chummage Cycler model trajectory 20 builds M according to double round Model B CM_LKinetic model.

Described multiple shooting method orbital exponent module 30 first aspect uses multiple shooting method respectively to M_BCMIt is modified, Obtain revised resonance Cycler dynamics of orbits model DM_BCM；Second aspect is to DM_BCMEmploying Fourth order Runge-Kutta obtains Resonance Cycler track SM_BCM；The third aspect uses multiple shooting method respectively to M_LIt is modified, obtains revised chummage Cycler dynamics of orbits model DM_L；Fourth aspect is to DM_LFourth order Runge-Kutta is used to obtain chummage Cycler track SM_L。

Described Lambert transfer orbit acquisition module 50 first aspect uses Gauss-global variable composite algorism to SM_BCM Process, obtain first set Lambert transfer orbitSecond aspect uses Differential correction algorithm pairIt is modified, obtains the revised Lambert of first set Transfer orbitThe third aspect is according to launch latitude and the white red angle of cut Changing Pattern is to SM_BCMCarry out launch window and the analysis of intersection window, obtain Spacecraft Launch and the opportunity entered the orbit；Four directions Face uses Gauss-global variable composite algorism to SM_LProcess, obtain the second set Lambert transfer orbit5th aspect uses differential correction algorithm pairIt is modified, after obtaining the second set correction Lambert transfer orbit6th aspect is by the Fuel Consumption of comparatively moon two-way mission, total Flight time can know SM_BCMAnd SM_LThe respective advantage of track, for the ground moon of survey of deep space of low-thrust spacecraft toward returning to one's post Business provides and optimizes design objective.

(1) resonance Cycler model trajectory 10 is built

In the present invention, obtain, according to circular re stricted three body problem MODEL C R3BP, the Cycler dynamics of orbits mould that resonates Type M_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 in figure 2 it is shown, ground moon barycenter inertial coodinate system O-XYZ and ground moon barycenter rotating coordinate system O-xyz.Wherein O is the ground moon Barycenter, with earth-moon system plane of movement for XY coordinate surface.In ground moon barycenter inertial coodinate system, X is that the initial time earth points to the moon Direction, Z-direction is the angular speed direction of earth-moon system, and Y-direction determines according to X and Z-direction right hand rule；Ground moon barycenter rotational coordinates Z direction in system with ground the moon barycenter inertial coodinate system Z-direction overlaps, x direction be always the earth sensing moon direction, y according to Right-handed system rule is set up.If the earth is P₁, the moon is P₂, spacecraft is P, and spacecraft P is under ground moon barycenter rotating coordinate system O-xyz Position be designated as (x_P,y_P,z_P)；Spacecraft P to earth P₁Distance be designated as Spacecraft P to moon P₂Distance be designated asSpacecraft P is to ground lunar geology The distance of heart O is designated as R.

Under ground moon barycenter rotating coordinate system O-xyz, the physical significance of formula (1) letter is:Represent on x direction Partial derivative；

U represents the potential function of spacecraft P, andμ₁Represent the moon and the earth Mass ratio, general value is 0.01215；

Represent the partial derivative on y direction；

Represent the partial derivative on z direction；

VA^xRepresent spacecraft P speed in the x direction；

VA^yRepresent spacecraft P speed in y-direction；

HA^xRepresent spacecraft P acceleration in the x direction；

HA^yRepresent spacecraft P acceleration in y-direction；

HA^zRepresent spacecraft P acceleration in a z-direction.

In the present invention, according to double round Model B CM to M_CR3BPIt is modified, obtains kinetic model under double circle Model B CM M_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)$

Under double round Model B CM, employ day ground rotating coordinate system O_S-X_SY_SZ_S, as shown in Figure 2 A, the sun is P₃, Spacecraft P is at day ground rotating coordinate system O_S-X_SY_SZ_SUnder position be designated as (xs_P,ys_P,zs_P)；Moon P₂Day rotate seat Mark system O_S-X_SY_SZ_SUnder position be designated asr_PSRepresent the nondimensionalization distance of spacecraft P and the sun, Andr_PERepresent the nondimensionalization distance of spacecraft P and the earth, andr_PMRepresent the nondimensionalization distance of spacecraft P and the moon, andμ₂Representing the mass ratio of the earth and the sun, general value is 0.000003003, m_MRepresent moon nondimensionalization quality.With day ground barycenter O_SFor initial point, day ground line be X_SAxle, points to the earth For X_SThe positive direction of axle.Y_SAxle is perpendicular to day ground plane of movement, according to the right-hand rule, it may be determined that Z_SThe direction of axle.

At day ground rotating coordinate system O_S-X_SY_SZ_SUnder, in formula (2), the physical significance of letter is:

Represent place M_BCMSpacecraft P under model is at X_SSpeed on direction；

Represent place M_BCMSpacecraft P under model is at Y_SSpeed on direction；

Represent place M_BCMSpacecraft P under model is at X_SAcceleration on direction；

Represent place M_BCMSpacecraft P under model is at Y_SAcceleration on direction；

Represent place M_BCMSpacecraft P under model is at Z_SAcceleration on direction.

(2) chummage Cycler model trajectory 20 is built

In the present invention, chummage Cycler dynamics of orbits model M is obtained according to double circle Model B CM_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)$

At day ground rotating coordinate system O_S-X_SY_SZ_SUnder, in formula (3), the physical significance of letter is:

Represent place M_LSpacecraft P under model is at X_SSpeed on direction；

Represent place M_LSpacecraft P under model is at Y_SSpeed on direction；

Represent place M_LSpacecraft P under model is at X_SAcceleration on direction；

Represent place M_LSpacecraft P under model is at Y_SAcceleration on direction；

Represent place M_LSpacecraft P under model is at Z_SAcceleration on direction.

(3) multiple shooting method orbital exponent module 30

In the present invention, use multiple shooting method to initialized M_BCMModel is modified, and obtains revised resonance Cycler dynamics of orbits model DM_BCM.Described DM_BCMGraphically characterize (as shown in Figure 3).Described multiple shooting method " calculating of a class flight optimization track " that reference is delivered for 1988, volume nine, the first phase, A21 to A22 page, author Wang Peide Deng.

M a cycle_BCMTaking the sampling point of 500 constant durations on model, any one sampling point is designated as ini_b (i) (i =1 ..., 500)；According to M_BCMModel uses Fourth order Runge-Kutta integration, obtains the position and speed quantity of state note of any point For ini_a (i) (i=1 ..., 500)；Ini_a (i :) is the matrix of 500 × 6.

In each time interval, choose 100 sampling points, and sampling point is used ODE45 integration, take out the position speed of end point Degree quantity of state φ (ini_a (i :)).The position and speed quantity of state φ (ini_a (i :) of note end point) with the position speed of any point The difference of degree quantity of state ini_a (i :) is F (i :), i.e. F (i :)=φ (ini_a (i :))-ini_a (i :).Application difference F (i :) it is so that resolving value and Practical Calculation value can be error free, should try one's best in simulation process and error is tapered to minimum.

Newton iteration method minimizing is used for above-mentioned difference F (i :).Described Newton iteration method was sent out with reference to 2011 " MATLAB of Newton iteration method realizes " of table, the 6th phase, page 20, author Yun Lei.It is 1 × 10 at setting accuracy^-10Condition Under, through 8 iteration, reach setting accuracy.

1st iteration	2.4925020681966167×10-1
		2nd iteration	9.5043845337684230×10-2
3rd iteration	2.3381313323453772×10-2
		The 4th iteration	1.9075592366893615×10-2
The 5th iteration	4.5915032506800547×10-4
		6th iteration	8.9114511750870658×10-5
7th iteration	4.5782899616269825×10-9
		8th iteration	8.9656802854146783×10-10

If more than set precision, continue to repeat above-mentioned iteration with Newton iteration method, until less than set precision Till.Thus, it is possible to obtain revised model DM_BCM。

To DM_BCMFourth order Runge-Kutta is used to obtain resonance Cycler track SM_BCM, this SM_BCMTrack is as shown in Figure 3.Figure In, each point has reacted spacecraft P position under ground moon barycenter inertial coodinate system O-XYZ in five cycles.Described quadravalence dragon lattice " runge kutta method and Mathematica thereof realize " that Ku Tafa reference is delivered for 2006, volume 18, the 2nd phase, the 72nd to the 73rd Page, author Chen Min.

In the present invention, use multiple shooting method to initialized M_LModel is modified, and obtains revised chummage Cycler dynamics of orbits model DM_L.Described DM_LGraphically characterize (as shown in Figure 3A).

M a cycle_LTaking the sampling point of 500 constant durations on model, any one sampling point is designated as ini_d (i) (i =1 ..., 500)；According to M_LModel uses Fourth order Runge-Kutta integration, and the position and speed quantity of state obtaining any point is designated as Ini_c (i) (i=1 ..., 500)；Ini_c (i :) is the matrix of 500 × 6.

In each time interval, choose 100 sampling points, and sampling point is used ODE45 integration, take out the position speed of end point Degree quantity of state θ (ini_c (i :)).The position and speed quantity of state θ (ini_c (i :) of note end point) with the position and speed of any point The difference of quantity of state ini_c (i :) is G (i :), i.e. G (i :)=θ (ini_c (i :))-ini_c (i :).Application difference G (i :) it is so that resolving value and Practical Calculation value can be error free, should try one's best in simulation process and error is tapered to minimum.

Newton iteration method minimizing is used for above-mentioned difference G (i :), is 1 × 10 at setting accuracy^-10Under the conditions of, Through 8 iteration, reach setting accuracy.

If more than set precision, continue to repeat above-mentioned iteration with Newton iteration method, until less than set precision Till.Thus, it is possible to obtain revised model DM_L。

To DM_LFourth order Runge-Kutta is used to obtain resonance Cycler track SM_L, this SM_LTrack is as shown in Figure 3A.In figure, Each point has reacted spacecraft P position under ground moon barycenter inertial coodinate system O-XYZ.

To DM_LFourth order Runge-Kutta is used to obtain chummage Cycler track SM_L。

(4) Lambert transfer orbit acquisition module 50

Lambert problem is the two_point boundary value problem in astrodynamics, in Spacecraft Rendezvous, missile intercept, interspace The fields such as navigation are widely used.As shown in Figure 4, origin endpoint E of spacecraft P₁, terminal E₂Position vector be respectively L₁And L₂, the focus of elliptical transfer orbit is positioned at the earth's core P₁.Origin endpoint E₁Time be designated as t₁, terminal E₂Time be designated as t₂, θ For angle of shift.

According to Lambert flight time theorem, on elliptical transfer orbit, the transfer time between any two points turns with oval Move the major semiaxis ra of track, half past two footpath sum (L₁+L₂) and central angle θ relevant, then it represents that for:

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

If L₁With L₂Sum is constant, and major semiaxis ra is constant, origin endpoint E₁With terminal E₂Between distance R_PFor constant, Then from origin endpoint E₁To terminal E₂Flight t transfer time_fAlso it is constant.The determination of elliptical transfer orbit and point-to-point transmission speed The selection of degree is the key of Lambert problem, and it is described as following Gauss problem: pursuit spacecraft P sets out and locates position arrow Amount (L₁) and velocity (v₁), spacecraft P terminal location vector (L₂) and velocity (v₂), flight transfer time is t_f, boat It device P is at elliptical transfer orbit E₁Initial velocity at Dian is v₁₀, spacecraft P is at elliptical transfer orbit E₂End speed at Dian For v₂₀.The target of this Gauss problem is to solve initial position speed increment Δ v₁With terminal location speed increment Δ v₂, enter And the impulse force size of the applying asked.

Gauss problem can be tried to achieve by following transcendental equations:

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

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

Lagrange coefficient one

Lagrange coefficient two

Lagrange coefficient three

Lagrange coefficient fourμ is time celestial body and primary body mass ratio, in round restriction In property problem of three bodies MODEL C R3BP, value is 0.01215；H is that spacecraft P is on the elliptical transfer orbit of various location Variable.

In formula (5), formula (6), it is known that L₁And v₁₀, or L₂And v₂₀It is assured that out that spacecraft P runs ellipse Circle transfer orbit.Obviously, once it is determined that Lagrange coefficient k, g,Gauss problem just can be readily solved.

In the present invention, global variable algorithm is used to solve for formula (5), formula (6).

" the Technique in Rendezvous and Docking mission planning " that described Gauss-global variable composite algorism reference is published for 2008, the 81st Page to page 100 content, page 142 to the 144th page, author Tang Guojin etc..

Assuming that the moon is 3 days to resonance Cycler inter-orbital transfer time, Lambert transfer orbit start point distance moonscape is high Degree 100km.Resonance Cycler track is 0.5 day to earth transfer time, and Lambert transfer orbit terminal is away from earth surface height 100km.Assuming that the earth is 0.5 day to resonance Cycler inter-orbital transfer time, Lambert transfer orbit start point distance earth surface is high Degree 100km.Resonance Cycler track is 1 day to moon transfer time, and Lambert transfer orbit terminal is away from moon earth surface height 100km。

Use Gauss-global variable composite algorism to SM_BCMProcess, obtain first set Lambert transfer orbitCorresponding Lambert track simulation result is as shown in Figure 5.In figure, in the ground moon 4 sections of dotted lines, Article 1 dotted line is had under barycenter inertial coodinate system O-XYZBy Cycler track to moon parking orbit Lambert transfer orbit；Article 2 dotted lineLambert transfer orbit by moon parking orbit to Cycler track；The Article three, dotted lineLambert transfer orbit by earth parking orbit to Cycler track；Article 4 dotted lineBy Cycler track is to the Lambert transfer orbit of earth parking orbit.Solid line is SM_BCMTrack.

Classical disome Lambert problem is typical two-point boundary value problem, uses the method for differential correction to be modified the A set of Lambert transfer orbitMakeover process is:

According to L₁、L₂Disome Lambert transfer orbit can be calculated with θ.Lambert becomes rail strategy turning to out position intersection Moving track, controlled quentity controlled variable is origin endpoint E₁The impulse speed increment of position, its deflection is designated as β；Differential correction algorithm will change Enter initial position speed increment Δ v₁And t transfer time_f, to realize end point E₂The intersection of position.Press close at the beginning of the iteration of true value Value, it is ensured that the convergence of differential correction algorithm iterative process.In order to spacecraft P is directed into target location L₂(i.e. theoretical whole End position vector), each iterative process is all by orbit integration to L₂Place, and t_fI.e. it is taken as this orbit integration time.Obviously, initially End points E₁Orbital velocity v of position₁₀Changes delta v₁Orbit integration time t will be caused_fChange, be designated as Δ t_f.Examine or check m Secondary iteration (is once designated as m+1 after the front m-1 that is once designated as that iteration is m time, iteration m time, then iteration is designated as m time when previous), right End point E₂The position vector of position is made single order Taylor and is launched, and can 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₁) represent that spacecraft P terminal location actual vector is about transfer time, at ellipse transfer rail Road E₁The function between initial velocity at Dian；

Represent terminal location vector and elliptical transfer orbit E₁The partial derivative between initial velocity at Dian, is reduced to

In the present invention, speed correction amount Δ v₁Should make:

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

Arrangement formula (7) and formula (8), 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)$

Represent the physical end position vector of spacecraft P.

Might as well takeAnd L₂There is identical abscissa component, at day ground rotating coordinate system O_S-X_SY_SZ_SUnder, and remember terminal Position vector and elliptical transfer orbit E₁The partial derivative between initial velocity at DianExpansion (9), can obtain:

$\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 represents day ground rotating coordinate system O_S-X_SY_SZ_SUnder Y_SThe position difference of axle；

δ z represents day ground rotating coordinate system O_S-X_SY_SZ_SUnder Z_SThe position difference of axle；

M₁₁Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian A line first row element；

M₁₂Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian A line secondary series element；

M₁₃Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian A line the 3rd column element；

M₂₁Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Two row first row elements；

M₂₂Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Two row secondary series elements；

M₂₃Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Two row the 3rd column elements；

M₃₁Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Three row first row elements；

M₃₂Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Three row secondary series elements；

M₃₃Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Three row the 3rd column elements；

Represent terminal location E₂Place terminates speed at X_SThe velocity component of axle；

Represent terminal location E₂Place terminates speed at Y_SThe velocity component of axle；

Represent terminal location E₂Place terminates speed at Z_SThe velocity component of axle；

Represent origin endpoint E₁Position initial velocity increment is at X_SThe velocity component of axle；

Represent origin endpoint E₁Position initial velocity increment is at Y_SThe velocity component of axle；

Represent origin endpoint E₁Position initial velocity increment is at Z_SThe velocity component of axle；

Δt_fRepresent origin endpoint E₁Orbital velocity v of position₁₀Increment Delta v₁The change of the orbit integration time caused Change.

OrderThen have:

$\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^TThe inversion of representing matrix C.

From the point of view of speed compared with a front m-1, then it is characterized as when the speed increment of previous m:

${\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)$

After setting accuracy, iterate, until meeting precision.Take off to resonance type Cycler for moon parking orbit The transfer orbit of track, setting accuracy is 1 × 10^-8Iteration result is as shown in the table.

1st iteration	2.0620622902499468×100
		2nd iteration	2.7641929155045797×10-1
3rd iteration	7.6981072705594636×10-2
		The 4th iteration	8.3862124467837824×10-3
The 5th iteration	1.7139690673584045×10-4
		6th iteration	1.0442627289698014×10-7

Through 6 iteration, reach setting accuracy.

For the transfer orbit of resonance type Cycler track to moon parking orbit, setting accuracy is 1 × 10^-8Iteration is tied The most as shown in the table.

1st iteration	2.5421929920011647×100
		2nd iteration	1.2546155043467997×10-1
3rd iteration	6.5637444323567223×10-2
		The 4th iteration	5.3832424435867889×10-3
The 5th iteration	2.2263784592810393×10-4
		6th iteration	2.0143221442658014×10-7

Through 6 iteration, reach setting accuracy.

Taking off to the transfer orbit of resonance type Cycler track for earth parking orbit, setting accuracy is 1 × 10^-8Repeatedly As shown in the table for result.

1st iteration	8.5792468327776916×10-1
		2nd iteration	7.9205098920078354×10-1
3rd iteration	1.4596634587845803×10-1
		The 4th iteration	2.5378649526028702×10-3
The 5th iteration	7.2114508959175548×10-7

Through 5 iteration, reach setting accuracy.

For the transfer orbit of resonance type Cycler track to earth parking orbit, setting accuracy is 1 × 10^-8Iteration is tied The most as shown in the table.

Through 7 iteration, reach setting accuracy.

In the present invention, differential correction algorithm pair is usedIt is modified, To first set revised Lambert transfer orbitSuch as Fig. 5 A, figure Shown in 5B: in figure, it is heavy line after correction.Spacecraft flight progress on actual transfer track, closer to heavy line, more accords with Resultant motion rule.Comparison diagram 5 and Fig. 5 A, whereinForCorrection,ForCorrection.Comparison diagram 5 with Fig. 5 B, whereinForCorrection,ForCorrection.

(1) analyze spacecraft to be taken off to the intersection of Cycler track by ground

First launch window analysis is carried out:

The most red (earth equatorial plane) angle of cut is in week monthly variation, the change during March 21 to 21 days March in 2006 in 2005 Scope is at [0 °, η], and wherein η value changes between 28.4 °～28.7 °.It is 27.9 ° that latitude is launched in Xichang, and latitude is launched in Wenchang Be 19 °, then there is launching opportunity space station can be emitted directly toward normal society face in both of which, and the window directly entered the orbit in Xichang is less than literary composition Prosperous；As shown in Figure 6, abscissa is the time, and ordinate is the most red (earth equatorial plane) angle of cut.For Wenchang, the width of launch window Degree is about 7 days, and interval is about 7 days；For Xichang, the width of launch window is about 2 days, and interval is about 13 days.In launch window In, suitable launching time (1 day 1 time) should be selected to be captured in normal society face by right ascension of ascending node.

Xichang launch latitude with reference within 2007, deliver " launching site, Xichang sketch China four of " firing dragon coming out of water " is big Satellite launch center (1) ", page 17, author Xiao Bo.

" the thunder and lightning environmental analysis of Wenchang, hainan rocket launching site " that the transmitting latitude reference of Wenchang is delivered for 2012, the 183rd Page, author is high.

Then intersection window analysis is carried out:

The Cycler track that resonates is transferred to, by parking orbit (height for 100km), in space station, take Huo Man branch mode or Lambert branch mode carries out the analysis of Fuel Consumption and transfer time, as shown in Figure 6 A and 6 B.For most saving combustion in Fig. 6 A The required speed increment of material Huo Man transfer is 3192m/s, and transfer time is 7 days.Fig. 6 B is the required speed of transfer of expense fuel Huo Man Degree increment is 4093.5m/s, and transfer time is 66.5min.

Owing to Fig. 6 A and Fig. 6 B reaches the apogean time roughly equal (about 7 days), and Fig. 6 A fuel is more saved, will As nominal intersection track: note nominal intersection be the 1st day launching track, if because nominal intersection launch be delayed, select the 2nd, 3, 4, within 5,6,7 days, launch, then required fuel is increased to 4093.5m/s by 3192km/s successively；Note nominal intersection is to send out for the 7th day Penetrate track, it is possible to select, from sky reciprocal number scale, within-6 ,-5 ,-4 ,-3 ,-2 ,-1 day, to launch, then required fuel successively by 3192km/s increases to 4093.5m/s；Note nominal intersection is the 1st～7 day launching track, then required fuel is successively by 3192km/s Increase to 4093.5m/s.The launching track of continuous 7 days as shown in Figure 6 C and 6 D shown in FIG., under ground moon barycenter inertial coodinate system O-XYZ As shown in Figure 6 C, under ground moon barycenter rotating coordinate system O-xyz, launching track is as shown in Figure 6 D for launching track.

(2) analyze spacecraft to be taken off to the intersection of resonance Cycler track by lunar surface

First launch window analysis is carried out:

The most red (moon equatorial plane) angle of cut is in week monthly variation, the change during March 21 to 21 days March in 2006 in 2005 Scope is at [0 °, ε], and wherein ε value changes between 6.64 °～6.86 °, as it is shown in fig. 7, abscissa is the time, ordinate is The most red (moon equatorial plane) angle of cut.Period of change is about 13.6 days.For the latitude landing point higher than 6.64 °, after lunar surface takes off Need to adjust orbit plane, resonance of being allowed for access Cycler track.As a example by south poles landing point, lunar surface takes off and enters 100km The ring moon SSO (Sun Synchronous Orbit)；Resonance Cycler orbital acquisition process is as follows: lifting apocynthion is to 20000km, and required speed increment is about 577m/ s；Apocynthion carries out inclination correction, and required speed increment is about 287.3m/s.The lifting of apocynthion can be next step resonance Cycler orbital rendezvous saves portion of energy.The moon revolution cycle is identical with the rotation period, only has every month twice chance can be by Right ascension of ascending node captures in normal society face.

Then intersection window analysis is carried out:

With the intersection time and resonance Cycler orbital phase as variable, using lunar surface parking orbit height 100km as constraint Condition, the search speed increment intersection window less than or equal to 1500m/s.Geographic logitude depends on right ascension of ascending node, according to " berthing Orbit altitude 100km " as constraints Search Results as shown in Figure 7 A, in fig. 7, T1 point, T2 point, T3 point and T4 point Represent the intersection situation meeting search condition and optimum.Geographic logitude as shown in Figure 7 B can be distributed in [-170 ° ,-135 °] ∪ [-43°,-5°]∪[45°,108°]∪[150°,165°].Landing point in above-mentioned interval, is expected to straight after being taken off by lunar surface Tap in normal society face；Landing point outside above-mentioned interval needs to wait that two weeks carries out right ascension of ascending node capture.Therefore, if selecting Specifying the landing point in region, lunar surface i.e. enters in normal society face after taking off: for resonance Cycler orbital phase [0.05,0.25] The interior intersection of can taking off at any time of ∪ [0.8,1].

In the present invention, use Gauss-global variable composite algorism to SM_LProcess, obtain the second set Lambert Transfer orbitTo SM_LProcessing method with to SM_BCMProcessing method be identical.In order to distinguish explanation, The relational expression of citation is added letter H distinguish.

According to Lambert flight time theorem, on elliptical transfer orbit, the transfer time between any two points turns with oval Move the major semiaxis Hra of track, half past two footpath sum (HL₁+HL₂) and central angle H θ relevant, then it represents that for:

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

If HL₁With HL₂Sum is constant, and major semiaxis Hra is constant, origin endpoint E₁With terminal E₂Between distance HR_PFor Constant, then from origin endpoint E₁To terminal E₂Flight Ht transfer time_fAlso it is constant.The determination of elliptical transfer orbit and two Between point, the selection of speed is the key of Lambert problem, and it is described as following Gauss problem: pursuit spacecraft P sets out place Position vector (HL₁) and velocity (Hv₁), spacecraft P terminal location vector (HL₂) and velocity (Hv₂), flight transfer Time is Ht_f, spacecraft P is at elliptical transfer orbit E₁Initial velocity at Dian is Hv₁₀, spacecraft P is at elliptical transfer orbit E₂Point The end speed at place is Hv₂₀.The target of this Gauss problem is to solve initial position speed increment Δ Hv₁With terminal location speed Degree increment Delta Hv₂, and then the impulse force size of the applying asked.

In the present invention, Gauss problem can be tried to achieve by following transcendental equations:

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

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

Lagrange coefficient one

Lagrange coefficient two

Lagrange coefficient three

Lagrange coefficient fourμ is time celestial body and primary body mass ratio, at double circles In Model B CM, value is 0.000003003；H is spacecraft P variable on the elliptical transfer orbit of various location.

In formula (14), formula (15), it is known that HL₁And Hv₁₀, or HL₂And Hv₂₀It is assured that out that spacecraft P transports The elliptical transfer orbit of row.Obviously, once it is determined that Lagrange coefficient k, g,Gauss problem just can be readily solved.

In the present invention, global variable algorithm is used to solve for formula (14), formula (15).

Assume that the earth to homoclinic orbit transfer time is 1 day, transfer orbit start point distance earth surface height 100km.Chummage Track to moon transfer time is 4 days, and transfer orbit terminal is away from moonscape height 100km.According to disome Lambert solution Method, corresponding chummage Lambert track simulation result is as seen in figure 7 c.

Use differential correction algorithm pairIt is modified, obtains the second set revised Lambert transfer TrackAfter second set correction, transfer orbit simulation result is as shown.

According to HL₁、HL₂Disome Lambert transfer orbit can be calculated with H θ.Lambert becomes rail strategy to out position intersection Transfer orbit, controlled quentity controlled variable is origin endpoint E₁The impulse speed increment of position, its deflection is designated as H β；Differential correction algorithm will Improve initial position speed increment Δ Hv₁And Ht transfer time_f, to realize end point E₂The intersection of position.Press close to changing of true value For initial value, it is ensured that the convergence of differential correction algorithm iterative process.In order to spacecraft P is directed into target location L₂(i.e. manage Opinion terminal location vector), each iterative process is all by orbit integration to HL₂Place, and Ht_fI.e. it is taken as this orbit integration time.Aobvious So, origin endpoint E₁Orbital velocity Hv of position₁₀Changes delta Hv₁Orbit integration time Ht will be caused_fChange, be designated as Δ Ht_f.(be once designated as m+1 after the front m-1 that is once designated as that iteration is m time, iteration m time, then iteration is designated as m time to examine or check the m time iteration When previous), to end point E₂The position vector of position is made single order Taylor and is launched, and can 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₁) represent that spacecraft P terminal location actual vector is about transfer time, at ellipse Transfer orbit E₁The function between initial velocity at Dian；

Represent terminal location vector and elliptical transfer orbit E₁The partial derivative between initial velocity at Dian, simplifies For

In the present invention, speed correction amount Δ Hv₁Should make:

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

Arrangement formula (16) and formula (17), can obtain:

${\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)$

Represent the physical end position vector of spacecraft P.

Might as well takeAnd L₂There is identical abscissa component, at day ground rotating coordinate system O_S-X_SY_SZ_SUnder, and remember terminal Position vector and elliptical transfer orbit E₁The partial derivative between initial velocity at DianExpansion (18), can :

$\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 day ground rotating coordinate system O_S-X_SY_SZ_SUnder Y_SThe position difference of axle；

δ Hz represents day ground rotating coordinate system O_S-X_SY_SZ_SUnder Z_SThe position difference of axle；

HM₁₁Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian A line first row element；

HM₁₂Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian A line secondary series element；

HM₁₃Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian A line the 3rd column element；

HM₂₁Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Two row first row elements；

HM₂₂Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Two row secondary series elements；

HM₂₃Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Two row the 3rd column elements；

HM₃₁Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Three row first row elements；

HM₃₂Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Three row secondary series elements；

HM₃₃Represent terminal location vector and elliptical transfer orbit E₁Of the partial derivative matrix between initial velocity at Dian Three row the 3rd column elements；

Represent terminal location E₂Place terminates speed at X_SThe velocity component of axle；

Represent terminal location E₂Place terminates speed at Y_SThe velocity component of axle；

Represent terminal location E₂Place terminates speed at Z_SThe velocity component of axle；

Represent origin endpoint E₁Position initial velocity increment is at X_SThe velocity component of axle；

Represent origin endpoint E₁Position initial velocity increment is at Y_SThe velocity component of axle；

Represent origin endpoint E₁Position initial velocity increment is at Z_SThe velocity component of axle；

ΔHt_fRepresent origin endpoint E₁Orbital velocity v of position₁₀Increment Delta v₁The change of the orbit integration time caused Change.

OrderThen have:

$\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^TThe inversion of representing matrix HC.

From the point of view of speed compared with a front m-1, then it is characterized as when the speed increment of previous m:

${\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)$

After setting accuracy, iterate, until meeting precision.Take off to resonance type Cycler for moon parking orbit The transfer orbit of track, setting accuracy is 1 × 10^-15Iteration result is as shown in the table.

1st iteration	2.8639968709823611×100
		2nd iteration	1.2827923060663518×100
3rd iteration	5.6581823976188161×10-1
		The 4th iteration	2.1973874406340671×10-1
The 5th iteration	1.5468617860521472×10-1
		6th iteration	4.9764011654783211×10-2
7th iteration	5.2642496344864611×10-3
		8th iteration	2.0846898885215647×10-5
9th iteration	2.2579833858870749×10-10
		10th iteration	6.5486562934038266×10-14

Through 10 iteration, reach setting accuracy.

For the transfer orbit of earth parking orbit to chummage type Cycler track, setting accuracy is 1 × 10^-11Iteration is tied The most as shown in the table.

1st iteration	4.0437979717451373×100
		2nd iteration	2.1369334069815822×100
3rd iteration	9.8469011235887069×10-2
		The 4th iteration	5.1247925851644960×10-4
The 5th iteration	5.5164314416937229×10-8
		6th iteration	3.3282881768155430×10-12

Through 6 iteration, reach setting accuracy.

In the present invention, compare based on resonance type Cycler track SM_BCMAnd chummage type Cycler track SM_LTurn by the ground moon The flight time (including becoming the rail time) of shifting task, fuel consumption.The general speed punching of rail is become with twice Lambert of each task Fuel Consumption levied by scale.

1) resonance type Cycler track scheme:

It it is 6.014 days by the total flight time of earth parking orbit to moon parking orbit.

It is 4.626km/s that twice Lambert during this period becomes the general speed momentum of rail.

2) chummage type Cycler track scheme:

It it is 20.695 days by the total flight time of earth parking orbit to moon parking orbit.

It is 4.450km/s that twice Lambert during this period becomes the general speed momentum of rail.

By contrast total flight time, total fuel consumption, it can be seen that the ground moon of based on Cycler track transfer scheme, Total flight time is longer, and total fuel consumption is less.This explanation, the ground moon of based on Cycler track the imagination of transfer scheme be feasible 's.Two class Cycler tracks have good periodicity.The cycle of resonance type Cycler track can set, and spacecraft can root Need according to task, in a planned way travel to and fro between the ground moon, the transmission carrying out goods and materials, personnel of rule.Chummage type Cycler track Cycle be also it is believed that and have passed through ground moon system LL₁Point, while completing ground moon turnaround mission, also may be used With moon system LL over the ground₁The scientific equipment arranged at Dian is replaced and overhauls.The ground moon based on this two classes Cycler track shifts The scheme flight time is longer, and total fuel consumption is less.Its best advantage of Appropriate application, can be designed that the ground moon comes and goes " public transport " system, for the survey of deep space of low-thrust spacecraft, the transmission of earth-moon system goods and materials has great meaning.

Claims (7)

1. the task simulation system come and gone by ground based on the Cycler track moon, it is characterised in that: this analogue system includes Build resonance Cycler model trajectory (10), build chummage Cycler model trajectory (20), multiple shooting method orbital exponent module And Lambert transfer orbit acquisition module (50) (30)；

Described resonance Cycler model trajectory (10) first aspect that builds builds according to circular re stricted three body problem MODEL C R3BP M_CR3BPKinetic model；Second aspect uses double circle Model B CM to described M_CR3BPKinetic model is optimized, and obtains M_BCM Kinetic model；

Described structure chummage Cycler model trajectory (20) builds M according to double circle Model B CM_LKinetic model；

Described multiple shooting method orbital exponent module (30) first aspect uses multiple shooting method respectively to M_BCMIt is modified, To revised resonance Cycler dynamics of orbits model DM_BCM；Second aspect is to DM_BCMFourth order Runge-Kutta is used to obtain altogether Shake Cycler track SM_BCM；The third aspect uses multiple shooting method respectively to M_LIt is modified, obtains revised chummage Cycler dynamics of orbits model DM_L；Fourth aspect is to DM_LFourth order Runge-Kutta is used to obtain chummage Cycler track SM_L；

Described Lambert transfer orbit acquisition module (50) first aspect uses Gauss-global variable composite algorism to SM_BCMEnter Row processes, and obtains first set Lambert transfer orbitSecond aspect uses micro- Divide correction algorithm pairIt is modified, obtains the revised Lambert of first set and turn Move trackThe third aspect becomes according to launch latitude and the white red angle of cut Law is to SM_BCMCarry out launch window and the analysis of intersection window, obtain Spacecraft Launch and the opportunity entered the orbit；Fourth aspect is adopted By Gauss-global variable composite algorism to SM_LProcess, obtain the second set Lambert transfer orbit 5th aspect uses differential correction algorithm pairIt is modified, obtains the second set revised Lambert transfer Track6th aspect by the Fuel Consumption of comparatively moon two-way mission, total flight time can Know SM_BCMAnd SM_LThe respective advantage of track, the ground moon two-way mission for the survey of deep space of low-thrust spacecraft provides optimization to set Meter index.

The task simulation system that ground based on the Cycler track moon the most according to claim 1 comes and goes, it is characterised in that: institute State in Lambert transfer orbit acquisition module (50), according to L₁、L₂Disome Lambert transfer orbit can be calculated with θ；Lambert Change rail strategy is to the transfer orbit of out position intersection, and controlled quentity controlled variable is origin endpoint E₁The impulse speed increment of position, its direction Angle is designated as β；Differential correction algorithm will improve initial position speed increment Δ v₁And t transfer time_f, to realize end point E₂Position The intersection at place；In order to spacecraft P is directed into target location L₂, each iterative process is all by orbit integration to L₂Place, and t_fI.e. take For this orbit integration time；Obviously, origin endpoint E₁Orbital velocity v of position₁₀Changes delta v₁When will cause orbit integration Between t_fChange, be designated as Δ t_f；Examine or check the m time iteration, to end point E₂The position vector of position is made single order Taylor and is launched, Can obtainRepresent spacecraft P Terminal location actual vector is about transfer time, at elliptical transfer orbit E₁The function between initial velocity at Dian；Table Show terminal location vector and elliptical transfer orbit E₁The partial derivative between initial velocity at Dian, is reduced toL₁For Origin endpoint E of spacecraft P₁Position vector；L₂Terminal E for spacecraft P₂Position vector；θ is angle of shift；

Speed correction amount Δ v₁L (t should be made_f+△t_f,v₁₀+△v₁)=L₂Set up, then have

Represent the physical end position vector of spacecraft P；

TakeAnd L₂There is identical abscissa component, at day ground rotating coordinate system O_S-X_SY_SZ_SUnder, and remember that terminal location vector is with ellipse Circle transfer orbit E₁The partial derivative between initial velocity at DianThen have

δ y represents day ground rotating coordinate system O_S-X_SY_SZ_SUnder Y_SThe position difference of axle；

δ z represents day ground rotating coordinate system O_S-X_SY_SZ_SUnder Z_SThe position difference of axle；

M₁₁Represent terminal location vector and elliptical transfer orbit E₁The first row of the partial derivative matrix between initial velocity at Dian First row element；

M₁₂Represent terminal location vector and elliptical transfer orbit E₁The first row of the partial derivative matrix between initial velocity at Dian Secondary series element；

M₁₃Represent terminal location vector and elliptical transfer orbit E₁The first row of the partial derivative matrix between initial velocity at Dian 3rd column element；

M₂₁Represent terminal location vector and elliptical transfer orbit E₁Second row of the partial derivative matrix between initial velocity at Dian First row element；

M₂₂Represent terminal location vector and elliptical transfer orbit E₁Second row of the partial derivative matrix between initial velocity at Dian Secondary series element；

M₂₃Represent terminal location vector and elliptical transfer orbit E₁Second row of the partial derivative matrix between initial velocity at Dian 3rd column element；

M₃₁Represent terminal location vector and elliptical transfer orbit E₁The third line of the partial derivative matrix between initial velocity at Dian First row element；

M₃₂Represent terminal location vector and elliptical transfer orbit E₁The third line of the partial derivative matrix between initial velocity at Dian Secondary series element；

M₃₃Represent terminal location vector and elliptical transfer orbit E₁The third line of the partial derivative matrix between initial velocity at Dian 3rd column element；

Represent terminal location E₂Place terminates speed at X_SThe velocity component of axle；

Represent terminal location E₂Place terminates speed at Y_SThe velocity component of axle；

Represent terminal location E₂Place terminates speed at Z_SThe velocity component of axle；

Represent origin endpoint E₁Position initial velocity increment is at X_SThe velocity component of axle；

Represent origin endpoint E₁Position initial velocity increment is at Y_SThe velocity component of axle；

Represent origin endpoint E₁Position initial velocity increment is at Z_SThe velocity component of axle；

△t_fRepresent origin endpoint E₁Orbital velocity v of position₁₀Increment △ v₁The change of the orbit integration time caused；

OrderThen haveC^TThe inversion of representing matrix C； From the point of view of speed compared with a front m-1, then it is characterized as when the speed increment of previous m

The task simulation system that ground based on the Cycler track moon the most according to claim 1 comes and goes, it is characterised in that: foundation Circular re stricted three body problem MODEL C R3BP obtains the Cycler dynamics of orbits model that resonates Wherein,

Represent the partial derivative on x direction；

U represents the potential function of spacecraft P, andμ₁Represent the matter of the moon and the earth Amount ratio, general value is 0.01215；R₁Represent spacecraft P to earth P₁Distance, R₂Represent spacecraft P to moon P₂Away from From；

Represent the partial derivative on y direction；

Represent the partial derivative on z direction；

VA^xRepresent spacecraft P speed in the x direction；

VA^yRepresent spacecraft P speed in y-direction；

HA^xRepresent spacecraft P acceleration in the x direction；

HA^yRepresent spacecraft P acceleration in y-direction；

HA^zRepresent spacecraft P acceleration in a z-direction.

The task simulation system that ground based on the Cycler track moon the most according to claim 1 comes and goes, it is characterised in that: depend on Chummage Cycler dynamics of orbits model is obtained according to double round Model B CMIts In,

Under double round Model B CM, employ day ground rotating coordinate system O_S-X_SY_SZ_S, the sun is P₃, spacecraft P is at day ground rotational coordinates It is O_S-X_SY_SZ_SUnder position be designated as (xs_P,ys_P,zs_P)；Moon P₂At day ground rotating coordinate system O_S-X_SY_SZ_SUnder position be designated asr_PSRepresent the nondimensionalization distance of spacecraft P and the sun, and r_PERepresent the nondimensionalization distance of spacecraft P and the earth, andr_PMRepresent Spacecraft P and the nondimensionalization distance of the moon, andμ₂Represent ground Ball and the mass ratio of the sun；m_MRepresent moon nondimensionalization quality；With day ground barycenter O_SFor initial point, day ground line be X_SAxle, points to The earth is X_SThe positive direction of axle；Y_SAxle is perpendicular to day ground plane of movement, according to the right-hand rule, determines Z_SThe direction of axle；

Represent place M_LSpacecraft P under model is at X_SSpeed on direction；

Represent place M_LSpacecraft P under model is at Y_SSpeed on direction；

Represent place M_LSpacecraft P under model is at X_SAcceleration on direction；

Represent place M_LSpacecraft P under model is at Y_SAcceleration on direction；

Represent place M_LSpacecraft P under model is at Z_SAcceleration on direction.

The task simulation system that ground based on the Cycler track moon the most according to claim 1 comes and goes, it is characterised in that: profit With launch latitude and white red angle of cut Changing Pattern gained, space station is transferred to, by parking orbit, the Cycler track that resonates, is taked Huo Man branch mode or Lambert branch mode carry out the analysis of Fuel Consumption and transfer time, and fuel saving Huo Man shifts Required speed increment is 3192m/s, and transfer time is 7 days；The required speed increment of transfer of expense fuel Huo Man is 4093.5m/s, Transfer time is 66.5min；Fuel saving Huo Man is shifted as nominal intersection track, if being delayed because nominal intersection is launched, Select within the 2nd, 3,4,5,6,7 days, launch, then required fuel is increased to 4093.5m/s by 3192km/s successively；Note nominal is handed over Can be the 7th day launching track, or select, from sky reciprocal number scale, within-6 ,-5 ,-4 ,-3 ,-2 ,-1 day, to launch, then required fuel 4093.5m/s is increased to successively by 3192km/s.

The task simulation system that ground based on the Cycler track moon the most according to claim 1 comes and goes, it is characterised in that: profit With launch latitude and white red angle of cut Changing Pattern gained, with intersection time and resonance Cycler orbital phase as variable, with the moon Face parking orbit height 100km is as constraints, the search speed increment intersection window less than or equal to 1500m/s；Geographical warp Degree depends on right ascension of ascending node, can search for the intersection feelings obtaining optimum as constraints according to parking orbit height 100km Condition.

The task simulation system that ground based on the Cycler track moon the most according to claim 1 comes and goes, it is characterised in that: ratio Relatively based on resonance type Cycler track SM_BCMAnd chummage type Cycler track SM_LThe ground moon transfer flight time of task, fuel Consume；The general speed momentum becoming rail with twice Lambert of each task characterizes Fuel Consumption；

Described resonance type Cycler track: be 6.014 days by the total flight time of earth parking orbit to moon parking orbit；? It is 4.626km/s that twice Lambert during this becomes the general speed momentum of rail；

Described chummage type Cycler track: be 20.695 days by the total flight time of earth parking orbit to moon parking orbit； It is 4.450km/s that twice Lambert during this period becomes the general speed momentum of rail.