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

Abstract

The invention discloses an earth-moon round trip task simulation system based on a Cycler orbit. The system is capable of simulating earth-moon round trip tasks, and the design of the system covers the design of spacecraft orbits, spacecraft task analysis, spacecraft dynamics simulation and the like. By establishing a resonant Cycler orbit model, a homoclinic Cycler orbit model, a multiple shooting algorithm orbit correction model and a Lambert transfer orbit acquisition model, the design of a resonant Cycler orbit and a homoclinic Cycler orbit is corrected, and the design of a Lambert transfer orbit is corrected, and the technical problems in launching rendezvous window analysis and earth-moon round strip system task analysis are solved. The system has the advantages that solar gravitational perturbation is considered when the earth-moon Cycler orbit is acquired, the approximate cyclic Cycler orbit is obtained by the multiple shooting algorithm; relative precision is relatively higher than that of the cyclic Cycler orbit computed by the prior methods.

Description

A kind of ground moon based on the Cycler track round task simulation system

Technical field

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

Background technology

The first Chinese moon exploration project is comprised of five big systemses such as moon probing satellite, carrier rocket, launching site, observing and controlling and Ground Application, and Chinese lunar exploration satellite engineering also has five large-engineering targets: the one, and first lunar exploration satellite of development and emission China; The 2nd, tentatively grasp the lunar orbiting exploration basic fundamental; The 3rd, carry out first lunar science and survey; The 4th, Primary Construction moon exploration space engineering system; The 5th, for the moon exploration successive projects is accumulated experience.The gordian technique of Gonna breakthrough moon probing satellite for this reason; Tentatively set up the large system of Chinese deep space exploration program; Every gordian techniquies such as checking useful load and data interpretation; Tentatively set up Chinese survey of deep space technology development system; Cultivate the corresponding talent team.

First Chinese moon exploration project four big science tasks are:

One, obtain the moonscape 3 D stereoscopic image, the essential structure of meticulous division moonscape and geomorphic unit, carry out the research of moonscape impact crater form, size, distribution, density etc., for division and the early stage history of evolution research at terrestrial planet surface age provides master data, and provide basic data etc. for lunar surface soft landing district's addressing and lunar base location optimization.

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

Three, survey lunar soil thickness, utilize microwave irradiation technology, obtain the thickness data of moonscape lunar soil, thereby obtain moonscape age and distribution thereof, and on this basis, content, resource distribution and the stock number etc. of estimation nuclear fusion fuel used to generate electricity helium 3.

Four, survey the space environment of the earth to the moon.The moon and earth mean distance are 380,000 kilometers, magnetic tail far away zone in the space, magnetic field of the earth, satellite, at this regional detectable solar cosmic ray high energy particle and solar wind plasma body, is studied the interaction of solar wind and the moon and magnetic field of the earth magnetic tail and the moon.

The Cycler track refers to and periodically travels to and fro between between the earth and the moon, near planet, is diversion and the track that do not stop.Running on aircraft on cyclic track can be for a long time (several years even more than ten years) keeps flight between planet and does not need orbit maneuver (or only needing very little orbit maneuver), thereby is considered to a kind of economical long-range mission mode of conserve energy.Cycler track scheme according to around day number of turns difference also can be subdivided into again whole circle track, half-turn track, general track etc.The 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, the CR3BP) motion of the relatively infinitesimal trisome of describing mass under the graviational interaction of two primary bodys that public barycenter moves in a circle around it.With reference in May, 2010, " Postgraduate School, Chinese Academy of Sciences's doctorate paper " of Li Mingtao, the related content of the 19th page to the 21st page.

Two round models (Bi-Circalar Model, BCM) are the basic models of the characteristics of motion of an infinitely small mass body of research under the universal gravitation effect of the system of the moon (the system for winding barycenter operates on circular orbit), the sun and the earth (operating on circular orbit around common barycenter).

Summary of the invention

The objective of the invention is to propose round task simulation system of a kind of ground moon based on the Cycler track, this system produces resonance type Cycler track and chummage type Cycler track under circular Restricted three-body model.Excentricity and the impact of other celestial bodies (as the sun, Jupiter etc.) gravitation due to lunar orbit, make the track of setting up under circular Restricted three-body model different from truth, and even, due to the impact of perturbation, it is constant that these tracks can't keep.Therefore, the track of setting up under circular Restricted three-body model of take is initial value, under two round models, is optimized, and utilizes 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 be regularly and the earth, the moon meet, produced and take the launch window that the earth or the moon is target.For completing complete Earth-moon transfer orbit, analogue system of the present invention is also according to the classical solution of disome Lambert problem, generated that the earth takes off, " earth-cycle transfer orbit " of moon landing; " the cycle transfer orbit-moon " that the moon takes off, the earth lands; And revise under the limbs model.Calculate and compare Fuel Consumption and the flight time of these two kinds of transfer orbits in analogue system of the present invention.

The present invention is round task simulation system of a kind of ground moon based on the Cycler track, and this analogue system includes and builds resonance Cycler model trajectory 10, builds chummage Cycler model trajectory 20, multiple shooting method track correcting module 30 and Lambert transfer orbit acquisition module 40.

Described structure resonance Cycler model trajectory 10 first aspects are according to CR3BP model construction M _cR3BPkinetic model; Second aspect adopts the BCM model to described M _cR3BPkinetic model is optimized, and obtains M _bCMkinetic model.

Described structure chummage Cycler model trajectory 20 is according to BCM model construction M _lkinetic model.

Described multiple shooting method track correcting module 30 first aspects adopt multiple shooting method respectively to M _bCMrevised, obtained revised resonance Cycler dynamics of orbits model DM _bCM; Second aspect is to DM _bCMadopt the fourth order Runge-Kutta method to obtain resonance Cycler track SM _bCM; The third aspect adopts multiple shooting method respectively to M _lrevised, obtained revised chummage Cycler dynamics of orbits model DM _l; Fourth aspect is to DM _ladopt the fourth order Runge-Kutta method to obtain chummage Cycler track SM _l.

Described Lambert transfer orbit acquisition module 40 first aspects adopt Gauss-global variable composite algorism to SM _bCMprocessed, obtained first set Lambert transfer orbit $S_{\mathrm{BCM}}=\{{\mathrm{<figure-callout\; id="1"\; label="s"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">s</figure-callout>}}_1^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="2"\; label="s"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">s</figure-callout>}}_2^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="3"\; label="s"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">s</figure-callout>}}_3^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="4"\; label="s"\; filenames="HDA0000382880340000041.png,HDA0000382880340000081.png"\; state="\{\{state\}\}">s</figure-callout>}}_4^{\mathrm{BCM}}\};$ Second aspect adopts differential correction algorithm pair revised, obtained the revised Lambert transfer orbit of first set ${\mathrm{DS}}_{\mathrm{BCM}}={\{\mathrm{<figure-callout\; id="1"\; label="ds"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">ds</figure-callout>}}_1^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="2"\; label="ds"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">ds</figure-callout>}}_2^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="3"\; label="ds"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">ds</figure-callout>}}_3^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="4"\; label="ds"\; filenames="HDA0000382880340000041.png,HDA0000382880340000081.png"\; state="\{\{state\}\}">ds</figure-callout>}}_4^{\mathrm{BCM}}\};$ The third aspect according to launch latitude and white red angle of cut Changing Pattern to SM _bCMcarry out the analysis of launch window and intersection window, the opportunity of obtaining Spacecraft Launch and entering the orbit; Fourth aspect adopts Gauss-global variable composite algorism to SM _lprocessed, obtained the second cover Lambert transfer orbit

the 5th aspect adopts differential correction algorithm pair

revised, obtained the revised Lambert transfer orbit of the second cover

the 6th aspect is by comparatively Fuel Consumption, the total flight time of month two-way mission can be known SM _bCMand SM _ltrack advantage separately, for the ground moon two-way mission of the survey of deep space of low-thrust spacecraft provides the optimal design 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.Not consider the impact of solar gravitation perturbation on Track desigh in order to revise in existing resonance type Cycler track and chummage type Cycle track.

2. system of the present invention has been considered the solar gravitation perturbation when acquisition ground moon Cycler track, and uses multiple shooting method to obtain approximate cycle Cycler track, compares precision with the cycle Cycler track of existing method calculating and increases.

3. system of the present invention is taken off while to the intersection window of Cycler track and spacecraft, being taken off intersection window to resonance Cycler track by lunar surface by ground analyzing spacecraft, a kind of iterative algorithm is provided, can, considering to obtain transition window in solar gravitation perturbation situation, compare precision with existing Lambert method result of calculation and increase.

The accompanying drawing explanation

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

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

Fig. 2 A is a day ground rotating coordinate system O _s-X _sy _sz _sthe coordinate system schematic diagram.

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

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

Fig. 4 is 2 schematic diagram that boundary value obtains in astrodynamics under the Lambert problem.

Fig. 5 is through Gauss of the present invention-global variable composite algorism treatment S M _bCMafter simulation result figure.

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

Fig. 5 B revises the simulation result figure of Lambert transfer orbit through revised the second cover of differential correction algorithm.

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

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

Fig. 6 B takes the track schematic diagram that fuel Huo Man shifts required speed increment most.

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

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

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

Fig. 7 A is the Search Results schematic diagram under the parking orbit height constraint condition that is 100km.

Fig. 7 B is the window schematic diagram of entering the orbit that geographic longitude is [170 ° ,-135 °] ∪ [43 ° ,-5 °] ∪ [45 °, 108 °] ∪ [150 °, 165 °].

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

Fig. 7 D is that the second cover after the differential correction algorithm is revised the Lambert transfer orbit ${\mathrm{DS}}_L={\{\mathrm{ds}}_1^L,{\mathrm{ds}}_2^L\}$ Simulation result.

Embodiment

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

Shown in Figure 1, the present invention is round task simulation system of a kind of ground moon based on the Cycler track, and this analogue system includes and builds resonance Cycler model trajectory 10, builds chummage Cycler model trajectory 20, multiple shooting method track correcting module 30 and Lambert transfer orbit acquisition module 40.

Described structure resonance Cycler model trajectory 10 first aspects are according to CR3BP model construction M _cR3BPkinetic model; Second aspect adopts the BCM model to described M _cR3BPkinetic model is optimized, and obtains M _bCMkinetic model.

Described structure chummage Cycler model trajectory 20 is according to BCM model construction M _lkinetic model.

Described multiple shooting method track correcting module 30 first aspects adopt multiple shooting method respectively to M _bCMrevised, obtained revised resonance Cycler dynamics of orbits model DM _bCM; Second aspect is to DM _bCMadopt the fourth order Runge-Kutta method to obtain resonance Cycler track SM _bCM; The third aspect adopts multiple shooting method respectively to M _lrevised, obtained revised chummage Cycler dynamics of orbits model DM _l; Fourth aspect is to DM _ladopt the fourth order Runge-Kutta method to obtain chummage Cycler track SM _l.

Described Lambert transfer orbit acquisition module 40 first aspects adopt Gauss-global variable composite algorism to SM _bCMprocessed, obtained first set Lambert transfer orbit $S_{\mathrm{BCM}}=\{{\mathrm{<figure-callout\; id="1"\; label="s"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">s</figure-callout>}}_1^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="2"\; label="s"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">s</figure-callout>}}_2^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="3"\; label="s"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">s</figure-callout>}}_3^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="4"\; label="s"\; filenames="HDA0000382880340000041.png,HDA0000382880340000081.png"\; state="\{\{state\}\}">s</figure-callout>}}_4^{\mathrm{BCM}}\};$ Second aspect adopts differential correction algorithm pair

revised, obtained the revised Lambert transfer orbit of first set ${\mathrm{DS}}_{\mathrm{BCM}}={\{\mathrm{<figure-callout\; id="1"\; label="ds"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">ds</figure-callout>}}_1^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="2"\; label="ds"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">ds</figure-callout>}}_2^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="3"\; label="ds"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">ds</figure-callout>}}_3^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="4"\; label="ds"\; filenames="HDA0000382880340000041.png,HDA0000382880340000081.png"\; state="\{\{state\}\}">ds</figure-callout>}}_4^{\mathrm{BCM}}\};$ The third aspect according to launch latitude and white red angle of cut Changing Pattern to SM _bCMcarry out the analysis of launch window and intersection window, the opportunity of obtaining Spacecraft Launch and entering the orbit; Fourth aspect adopts Gauss-global variable composite algorism to SM _lprocessed, obtained the second cover Lambert transfer orbit

the 5th aspect adopts differential correction algorithm pair

revised, obtained the revised Lambert transfer orbit of the second cover

the 6th aspect is by comparatively Fuel Consumption, the total flight time of month two-way mission can be known SM _bCMand SM _ltrack advantage separately, for the ground moon two-way mission of the survey of deep space of low-thrust spacecraft provides the optimal design index.

(1) build resonance Cycler model trajectory 10

In the present invention, according to the CR3BP model Cycler dynamics of orbits model M that obtains resonating _cR3BP:

$M_{\mathrm{CR}3\mathrm{BP}}=\left\{\begin{array}{c}\frac{\&PartialD;U}{\&PartialD;x}={\mathrm{HA}}^x-2V{A}^y\\ \frac{\&PartialD;U}{\&PartialD;y}={\mathrm{HA}}^y+2V{A}^x\\ \frac{\&PartialD;U}{\&PartialD;z}={\mathrm{HA}}^z\end{array}\right.---\left(1\right)$

As shown in Figure 2, ground moon barycenter inertial coordinates system O-XYZ and ground month barycenter rotating coordinate system O-xyz.Wherein O is ground month barycenter, take the earth-moon system plane of movement as the XY coordinate surface.In ground month barycenter inertial coordinates system, X is the direction that the initial time earth points to the moon, the angular velocity direction that the Z direction is earth-moon system, and Y-direction is determined according to X and Z direction right hand rule; Z direction in ground month barycenter rotating coordinate system overlaps with the Z direction of ground month barycenter inertial coordinates system, and the x direction is always the direction that the earth points to the moon, and y sets up according to the right-handed system rule.If the earth is P ₁, the moon is P ₂, spacecraft is P, the position of spacecraft P under ground moon barycenter rotating coordinate system O-xyz is designated as (x _p, y _p, z _p); Spacecraft P is to earth P ₁distance be designated as R ₁( $R_1={\&lsqb;{(x_P-{\mathrm{\&mu;}}_1)}^2+{y_P}^2+{z_P}^2\&rsqb;}^{\frac{1}{2}}$ ), spacecraft P is to moon P ₂distance be designated as R ₂( $R_2={\&lsqb;{(x_P+1-{\mathrm{\&mu;}}_1)}^2+{y_P}^2+{z_{\mathrm{<figure-callout\; id="2"\; label="P"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">P</figure-callout>}}}^2\&rsqb;}^{\frac{1}{2}}$ ), the distance of spacecraft P month barycenter O to ground is designated as R.

Under ground moon barycenter rotating coordinate system O-xyz, the physical significance of formula (1) letter is:

mean the partial derivative on the x direction;

U means the potential function of spacecraft P, and $U=\frac{1}{2}({x_{\mathrm{<figure-callout\; id="2"\; label="P"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">P</figure-callout>}}}^2+{x_P}^2)+\frac{1-{\mathrm{\&mu;}}_1}{{\mathrm{<figure-callout\; id="1"\; label="R"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">R</figure-callout>}}_1}+\frac{{\mathrm{\&mu;}}_1}{{\mathrm{<figure-callout\; id="2"\; label="R"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">R</figure-callout>}}_2},$ μ ₁the mass ratio that means the moon and the earth, general value is 0.01215;

mean the partial derivative on the y direction;

mean the partial derivative on the z direction;

VA ^xmean the speed of spacecraft P on the x direction;

VA ^ymean the speed of spacecraft P on the y direction;

HA ^xmean the acceleration of spacecraft P on the x direction;

HA ^ymean the acceleration of spacecraft P on the y direction;

HA ^zmean the acceleration of spacecraft P on the z direction.

In the present invention, according to the BCM model to M _cR3BPrevised, obtained kinetic model M under the BCM model _bCM:

$M_{\mathrm{BCM}}=\left\{\begin{array}{c}{\mathrm{HB}}^{X_S}={2\mathrm{VB}}^{Y_S}+{\mathrm{xs}}_P-(1-{\mathrm{\&mu;}}_2)\frac{{\mathrm{xs}}_P+{\mathrm{\&mu;}}_2}{{\mathrm{<figure-callout\; id="3"\; label="r\; PS"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{PS}}^{}}\\ {\mathrm{\&mu;}}_2\frac{{\mathrm{xs}}_P-1+{\mathrm{\&mu;}}_2}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{xs}}_P-x_{P_2}}{{\mathrm{<figure-callout\; id="3"\; label="r\; PM"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{PM}}^{}}\\ {\mathrm{HB}}^{Y_S}=-2V{B}^{X_S}+{\mathrm{ys}}_P-(1-{\mathrm{\&mu;}}_2)\frac{{\mathrm{ys}}_P}{r_{\mathrm{PS}}^3}-\\ {\mathrm{\&mu;}}_2\frac{{\mathrm{ys}}_P}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{ys}}_P-y_{P_2}}{{\mathrm{<figure-callout\; id="3"\; label="r\; PM"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{PM}}^{}}\\ {\mathrm{HB}}^{Z_S}=-(1-{\mathrm{\&mu;}}_2)\frac{{\mathrm{zs}}_P}{r_{\mathrm{PS}}^3}-{\mathrm{\&mu;}}_2\frac{{\mathrm{zs}}_P}{r_{\mathrm{PM}}^3}-m_M\frac{{\mathrm{zs}}_P}{r_{\mathrm{PM}}^3}\end{array}\right.---\left(2\right)$

Under the BCM model, used a day ground rotating coordinate system P _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 ₂at day ground rotating coordinate system O _s-X _sy _sz _sunder position be designated as r _pSthe nondimensionalization distance that means spacecraft P and the sun, and $r_{\mathrm{PS}}={[{({\mathrm{xs}}_P+{\mathrm{\&mu;}}_2)}^2+{{\mathrm{ys}}_P}^2+{{\mathrm{zs}}_P}^2]}^{\frac{1}{2}};$ R _pEthe nondimensionalization distance that means spacecraft P and the earth, and $r_{\mathrm{PE}}={[{({\mathrm{xs}}_P-1+{\mathrm{\&mu;}}_2)}^2+{{\mathrm{ys}}_P}^2+{{\mathrm{zs}}_P}^2]}^{\frac{1}{2}};$ R _pMthe nondimensionalization distance that means spacecraft P and the moon, and $r_{\mathrm{PM}}={[{({\mathrm{xs}}_P-x_{P_2})}^2+{({\mathrm{ys}}_P-y_{P_2})}^2+{{\mathrm{zs}}_P}^2]}^{\frac{1}{2}};$ μ ₂the mass ratio that means the earth and the sun, general value is 0.000003003, m _mmean moon nondimensionalization quality.With day ground barycenter O _sfor initial point, a day ground line is X _saxle, the sensing earth is X _sthe positive dirction of axle.Y _saxle, perpendicular to day ground plane of movement, according to the right-hand rule, can be determined 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:

mean place M _bCMspacecraft P under model is at X _sspeed on direction;

mean place M _bCMspacecraft P under model is at Y _sspeed on direction;

mean place M _bCMspacecraft P under model is at X _sacceleration on direction;

mean place M _bCMspacecraft P under model is at Y _sacceleration on direction;

mean place M _bCMspacecraft P under model is at Z _sacceleration on direction.

(2) build chummage Cycler model trajectory 20

In the present invention, obtain chummage Cycler dynamics of orbits model M according to the BCM model _l:

$M_{\mathrm{BCM}}=\left\{\begin{array}{c}{\mathrm{HB}}^{X_S}={2\mathrm{VB}}^{Y_S}+{\mathrm{xs}}_P-(1-{\mathrm{\&mu;}}_2)\frac{{\mathrm{xs}}_P+{\mathrm{\&mu;}}_2}{{\mathrm{<figure-callout\; id="3"\; label="r\; PS"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{PS}}^{}}\\ {\mathrm{\&mu;}}_2\frac{{\mathrm{xs}}_P-1+{\mathrm{\&mu;}}_2}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{xs}}_P-x_{P_2}}{{\mathrm{<figure-callout\; id="3"\; label="r\; PM"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{PM}}^{}}\\ {\mathrm{HB}}^{Y_S}=-2V{B}^{X_S}+{\mathrm{ys}}_P-(1-{\mathrm{\&mu;}}_2)\frac{{\mathrm{ys}}_P}{r_{\mathrm{PS}}^3}-\\ {\mathrm{\&mu;}}_2\frac{{\mathrm{ys}}_P}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{ys}}_P-y_{P_2}}{{\mathrm{<figure-callout\; id="3"\; label="r\; PM"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">r</figure-callout>}}_{\mathrm{PM}}^{}}\\ {\mathrm{HB}}^{Z_S}=-(1-{\mathrm{\&mu;}}_2)\frac{{\mathrm{zs}}_P}{r_{\mathrm{PS}}^3}-{\mathrm{\&mu;}}_2\frac{{\mathrm{zs}}_P}{r_{\mathrm{PM}}^3}-m_M\frac{{\mathrm{zs}}_P}{r_{\mathrm{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:

mean place M _lspacecraft P under model is at X _sspeed on direction;

mean place M _lspacecraft P under model is at Y _sspeed on direction;

mean place M _lspacecraft P under model is at X _sacceleration on direction;

mean place M _lspacecraft P under model is at Y _sacceleration on direction;

mean place M _lspacecraft P under model is at Z _sacceleration on direction.

(3) multiple shooting method track correcting module 30

In the present invention, adopt multiple shooting method to initialized M _bCMmodel is revised, and obtains revised resonance Cycler dynamics of orbits model DM _bCM.Described DM _bCMwith figure (as shown in Figure 3), characterized.Described multiple shooting method is with reference to " calculating of a class flight optimization track " within 1988, delivered, the 9th volume, the first phase, A21 to A22 page, author Wang Peide etc.

M at one-period _bCMget the sampling point of 500 constant durations on model, any sampling point be designated as ini_b (i) (i=1 ..., 500); According to M _bCMmodel adopts fourth order Runge-Kutta method integration, the position and speed quantity of state that obtains any point be designated as ini_a (i) (i=1 ..., 500); Ini_a (i :) is the matrix of 500 * 6.

On each time interval, choose 100 sampling points, and sampling point is adopted to the ODE45 integration, take out the position and speed quantity of state of end point

the position and speed quantity of state of note end point

with the difference of the position and speed quantity of state ini_a of any point (i :) be F (i :),

application difference F (i :) makes the value of resolving and the actual computation value can be error free, in simulation process, should error be tapered to minimum as far as possible.

Adopt Newton iteration method minimizing for above-mentioned difference F (i :).Described Newton iteration method is with reference to " MATLAB of Newton iteration method realizes " of within 2011, delivering, the 6th phase, the 20th page, author Yun Lei.At setting accuracy, be 1 * 10 ^-10under condition, through 8 iteration, reach setting accuracy.

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

The precision set if be greater than, continue to repeat above-mentioned iteration with Newton iteration method, until be less than the precision set.Thus, can obtain revised model DM _bCM.

To DM _bCMadopt the fourth order Runge-Kutta method to obtain resonance Cycler track SM _bCM, this SM _bCMtrack as shown in Figure 3.In figure, each point has reacted the position of spacecraft P under ground moon barycenter inertial coordinates system O-XYZ in five cycles." runge kutta method and Mathematica thereof realize " that described fourth order Runge-Kutta method reference is delivered in 2006, the 18th volume, the 2nd phase, the 72nd to the 73rd page, author Chen Chi is quick.

In the present invention, adopt multiple shooting method to initialized M _lmodel correction obtains revised chummage Cycler dynamics of orbits model DM _l.Described DM _lwith figure (as shown in Figure 3A), characterized.

M at one-period _lget the sampling point of 500 constant durations on model, any sampling point be designated as ini_d (i) (i=1 ..., 500); According to M _lmodel adopts fourth order Runge-Kutta method integration, the position and speed quantity of state that obtains any point be designated as ini_c (i) (i=1 ..., 500); Ini_c (i :) is the matrix of 500 * 6.

On each time interval, choose 100 sampling points, and sampling point is adopted to the ODE45 integration, the position and speed quantity of state θ of taking-up end point (ini_c (i :)).The position and speed quantity of state θ of note end point (ini_c (i :)) with the difference of the position and speed quantity of state ini_c of any point (i :), be G (i, :), i.e. G (i :)=θ (ini_c (i, :))-ini_c (i :).Application difference G (i :) makes the value of resolving and the actual computation value can be error free, in simulation process, should error be tapered to minimum as far as possible.

Adopting Newton iteration method minimizing for above-mentioned difference G (i :), is 1 * 10 at setting accuracy ^-10under condition, through 8 iteration, reach setting accuracy.

The 1st iteration	1.4925020681966167×10 -1
		The 2nd iteration	8.5043845337684230×10 -2
The 3rd iteration	2.3461313323453772×10 -2
		The 4th iteration	1.9075592366893615×10 -3

The 5th iteration	5.5915036547800547×10 -4
		The 6th iteration	8.9114511750870658×10 -7
The 7th iteration	4.5782899616269825×10 -9
		The 8th iteration	5.1435802854146783×10 -10

The precision set if be greater than, continue to repeat above-mentioned iteration with Newton iteration method, until be less than the precision set.Thus, can obtain revised model DM _l.

To DM _ladopt the fourth order Runge-Kutta method to obtain resonance Cycler track SM _l, this SM _ltrack as shown in Figure 3A.In figure, each point has reacted the position of spacecraft P under ground moon barycenter inertial coordinates system O-XYZ.

To DM _ladopt the fourth order Runge-Kutta method to obtain chummage Cycler track SM _l.

(4) Lambert transfer orbit acquisition module 40

The Lambert problem is 2 boundary value problems in astrodynamics, in the fields such as Spacecraft Rendezvous, missile intercept, interplanetary flight, is widely used.As shown in Figure 4, the 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 ₂, θ is angle of shift.

According to Lambert flight time theorem, the major semi-axis ra of the transfer time on elliptical transfer orbit between any two points and elliptical transfer orbit, half past two footpath sum (L ₁+ L ₂) and central angle θ relevant, be expressed as:

t _f=W(ra，(L ₁+L ₂)，R _P) (4)

If L ₁with L ₂sum is constant, and major semi-axis ra is constant, origin endpoint E ₁with terminal E ₂between distance R _pfor constant, from origin endpoint E ₁to terminal E ₂flight t transfer time _fit is also constant.Definite and selection poimt-to-point speed of elliptical transfer orbit is the key of Lambert problem, and it is described to following Gauss's problem: pursuit spacecraft P sets out and locates position vector (L ₁) and velocity (v ₁), spacecraft P terminal location vector (L ₂) and velocity (v ₂), flight transfer time is t _f, spacecraft P is at elliptical transfer orbit E ₁the initial velocity at some place is v ₁₀, spacecraft P is at elliptical transfer orbit E ₂the end speed at some place is v ₂₀.The target of this Gauss's problem is in order to solve initial position speed increment Δ v ₁with terminal location speed increment Δ v ₂, and then the impulse force size applied of asking.

For Gauss's problem, can be tried to achieve by following transcendental equation group:

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

${\mathrm{<figure-callout\; id="20"\; label="v"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">v</figure-callout>}}_{20}=\stackrel{\&CenterDot;}{k}\&times;{\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">L</figure-callout>}}_1+\stackrel{\&CenterDot;}{g}\&times;v_{10}---\left(6\right)$

Lagrange coefficient one $k=1-\frac{{\mathrm{\&mu;L}}_2}{h^2}(1-\cos \mathrm{\&theta;});$

Lagrange coefficient two $\stackrel{\&CenterDot;}{k}=\frac{\mathrm{\&mu;}}{h}\frac{1-\cos \mathrm{\&theta;}}{\sin \mathrm{\&theta;}}\&lsqb;\frac{\mathrm{\&mu;}}{h}(1-\cos \mathrm{\&theta;})-\frac{1}{L_1}-\frac{1}{{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}\&rsqb;;$

Lagrange coefficient three

Lagrange coefficient four

μ is time celestial body and primary body mass ratio, and in the CR3BP model, value is 0.01215; H is the variable of spacecraft P on the elliptical transfer orbit at diverse location place.

In formula (5), formula (6), known L ₁and v ₁₀, or L ₂and v ₂₀just can determine the elliptical transfer orbit of spacecraft P operation.Obviously, once determined Lagrangian coefficient k, g, gauss's problem just can be readily solved.

In the present invention, for formula (5), formula (6), adopt the global variable algorithm to be solved.

Described Gauss-global variable composite algorism is with reference to " the Technique in Rendezvous and Docking mission planning " of publication in 2008, the 81st page to the 100th page content, the 142nd page to the 144th page, author Tang Guojin etc.

Suppose that the moon is 3 days to the Cycler inter-orbital transfer time of resonating, Lambert transfer orbit start point distance moonscape height 100km.Be 0.5 day resonance Cycler track to earth transfer time, and Lambert transfer orbit terminal is apart from earth surface height 100km.Suppose that the earth is 0.5 day to the Cycler inter-orbital transfer time of resonating, Lambert transfer orbit start point distance earth surface height 100km.Be 1 day resonance Cycler track to moon transfer time, and Lambert transfer orbit terminal is apart from moon earth surface height 100km.

Adopt Gauss-global variable composite algorism to SM _bCMprocessed, obtained first set Lambert transfer orbit $S_{\mathrm{BCM}}=\{{\mathrm{<figure-callout\; id="1"\; label="s"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">s</figure-callout>}}_1^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="2"\; label="s"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">s</figure-callout>}}_2^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="3"\; label="s"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">s</figure-callout>}}_3^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="4"\; label="s"\; filenames="HDA0000382880340000041.png,HDA0000382880340000081.png"\; state="\{\{state\}\}">s</figure-callout>}}_4^{\mathrm{BCM}}\},$ Corresponding Lambert track simulation result as shown in Figure 5.In figure, under ground moon barycenter inertial coordinates system O-XYZ, 4 sections dotted lines are arranged, article one dotted line lambert transfer orbit by the Cycler track to moon parking orbit; The second dotted line

lambert transfer orbit by moon parking orbit to the Cycler track; Article three, dotted line

lambert transfer orbit by earth parking orbit to the Cycler track; Article four, dotted line lambert transfer orbit by the Cycler track to earth parking orbit.Solid line is SM _bCMtrack.

Classical disome Lambert problem is typical two-point boundary value problem, adopts the method for differential correction to be revised first set Lambert transfer orbit

makeover process is:

According to L ₁, L ₂can calculate disome Lambert transfer orbit with θ.Lambert becomes the transfer orbit of rail strategy to the out position intersection, and controlled quentity controlled variable is origin endpoint E ₁the impulse speed increment of position, its deflection is designated as β; The differential correction algorithm will improve initial position speed increment Δ v ₁reach t transfer time _f, to realize end point E ₂the intersection of position.Press close to the iterative initial value of true value, can guarantee the convergence of differential correction algorithm iterative process.For spacecraft P is directed into to target location L ₂(being theoretical terminal location vector), each iterative process all by orbit integration to L ₂place, and t _fbe taken as this orbit integration time.Obviously, origin endpoint E ₁the orbital velocity v of position ₁₀changes delta v ₁to cause orbit integration time t _fvariation, be designated as Δ t _f.Examine or check the m time iteration (the front m-1 that once is designated as that iteration is m time, once be designated as m+1 after iteration m time, iteration be designated as for m time ought be last time), to end point E ₂the position vector of position is made single order Taylor and is launched, and can obtain:

$L(t_f+\mathrm{\&Delta;}{t}_f,{\mathrm{<figure-callout\; id="10"\; label="v"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">v</figure-callout>}}_{10}+\mathrm{\&Delta;}{v}_1)={\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">L</figure-callout>}}_2+\frac{{\&PartialD;\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}{{\&PartialD;\mathrm{<figure-callout\; id="10"\; label="v"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">v</figure-callout>}}_{10}}\&times;\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">v</figure-callout>}}_1+{\mathrm{<figure-callout\; id="20"\; label="v"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">v</figure-callout>}}_{20}\&times;\mathrm{\&Delta;}{t}_f---\left(7\right)$

L(t _f+ Δ t _f, v ₁₀+ Δ v ₁) mean that spacecraft P terminal location actual vector is about transfer time, at elliptical transfer orbit E ₁function between the initial velocity at some place;

mean terminal location vector and elliptical transfer orbit E ₁partial derivative between the initial velocity at some place, be reduced to $M=\frac{{\&PartialD;\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}{{\&PartialD;\mathrm{<figure-callout\; id="10"\; label="v"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">v</figure-callout>}}_{10}};$

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{\&delta;}{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">L</figure-callout>}}_2=L_2-\stackrel{\&OverBar;}{{\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}=\frac{{\&PartialD;\mathrm{<figure-callout\; id="2"\; label="L"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">L</figure-callout>}}_2}{{\&PartialD;\mathrm{<figure-callout\; id="10"\; label="v"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">v</figure-callout>}}_{10}}\&times;\mathrm{\&Delta;}{v}_1-{\mathrm{<figure-callout\; id="20"\; label="v"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">v</figure-callout>}}_{20}\&times;\mathrm{\&Delta;}{t}_f---\left(9\right)$

the physical end position vector that means spacecraft P.

Might as well get

and L ₂there is identical horizontal ordinate component, at day ground rotating coordinate system O _s-X _sy _sz _sunder, and note terminal location vector and elliptical transfer orbit E ₁partial derivative between the initial velocity at some place

expansion (9) can obtain:

$\left[\begin{array}{c}0\\ \mathrm{\&delta;y}\\ \mathrm{\&delta;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]\&times;\left[\begin{array}{c}\mathrm{\&Delta;}{v}_1^x\\ \mathrm{\&Delta;}{v}_1^y\\ {\mathrm{\&Delta;v}}_1^z\\ {\mathrm{\&Delta;t}}_f\end{array}\right]---\left(10\right)$

δ y means a day ground rotating coordinate system O _s-X _sy _sz _sunder Y _sthe position difference of axle;

δ z means a day ground rotating coordinate system O _s-X _sy _sz _sunder Z _sthe position difference of axle;

M ₁₁mean terminal location vector and elliptical transfer orbit E ₁the first row first row element of the partial derivative matrix between the initial velocity at some place;

M ₁₂mean terminal location vector and elliptical transfer orbit E ₁the first row secondary series element of the partial derivative matrix between the initial velocity at some place;

M ₁₃mean terminal location vector and elliptical transfer orbit E ₁the first row the 3rd column element of the partial derivative matrix between the initial velocity at some place;

M ₂₁mean terminal location vector and elliptical transfer orbit E ₁the second row first row element of the partial derivative matrix between the initial velocity at some place;

M ₂₂mean terminal location vector and elliptical transfer orbit E ₁the second row secondary series element of the partial derivative matrix between the initial velocity at some place;

M ₂₃mean terminal location vector and elliptical transfer orbit E ₁the second row the 3rd column element of the partial derivative matrix between the initial velocity at some place;

M ₃₁mean terminal location vector and elliptical transfer orbit E ₁the third line first row element of the partial derivative matrix between the initial velocity at some place;

M ₃₂mean terminal location vector and elliptical transfer orbit E ₁the third line secondary series element of the partial derivative matrix between the initial velocity at some place;

M ₃₃mean terminal location vector and elliptical transfer orbit E ₁the third line the 3rd column element of the partial derivative matrix between the initial velocity at some place;

mean terminal location E ₂end speed in place's is at X _sthe speed component of axle;

mean terminal location E ₂end speed in place's is at Y _sthe speed component of axle;

mean terminal location E ₂end speed in place's is at Z _sthe speed component of axle;

mean origin endpoint E ₁position initial velocity increment is at X _sthe speed component of axle;

mean origin endpoint E ₁position initial velocity increment is at Y _sthe speed component of axle;

mean origin endpoint E ₁position initial velocity increment is at Z _sthe speed component of axle;

Δ t _fmean origin endpoint E ₁the orbital velocity v of position ₁₀increment Delta v ₁the variation of the orbit integration time caused.

Order $C=\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],$ Have:

$\left[\begin{array}{c}{\mathrm{\&Delta;v}}_1^x\\ {\mathrm{\&Delta;v}}_1^y\\ {\mathrm{\&Delta;v}}_1^z\\ {\mathrm{\&Delta;t}}_f\end{array}\right]=-{\left(C^{T}C\right)}^{-1}{\mathrm{<figure-callout\; id="0"\; label="C\; T"\; filenames="HDA0000382880340000012.png,HDA0000382880340000022.png"\; state="\{\{state\}\}">C</figure-callout>}}^T\left[\begin{array}{c}0\\ \mathrm{\&delta;y}\\ \mathrm{\&delta;z}\end{array}\right]---\left(11\right)$

C ^tthe inversion of representing matrix C.

Therefore, with respect to the speed of a front m-1, last time the speed increment of m is characterized by:

$\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">v</figure-callout>}}_1=\left[\begin{array}{ccc}{\mathrm{\&Delta;v}}_1^x& {\mathrm{\&Delta;v}}_1^y& {\mathrm{\&Delta;v}}_1^z\end{array}\right]---\left(12\right)$

After setting accuracy, iterate, until meet precision.For moon parking orbit, take off to the transfer orbit of resonance type Cycler track, setting accuracy is 1 * 10 ^-8the iteration result is as shown in the table.

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

Through 6 iteration, reach setting accuracy.

Transfer orbit for resonance type Cycler track to moon parking orbit, setting accuracy is 1 * 10 ^-8the iteration result is as shown in the table.

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

Through 6 iteration, reach setting accuracy.

For earth parking orbit, take off to the transfer orbit of resonance type Cycler track, setting accuracy is 1 * 10 ^-8the iteration result is as shown in the table.

The 1st iteration	8.5792468327776916×10 -1
		The 2nd iteration	7.9205098920078354×10 -1
The 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.

Transfer orbit for resonance type Cycler track to earth parking orbit, setting accuracy is 1 * 10 ^-8the iteration result is as shown in the table.

The 1st iteration	8.2182773945012833×10 -1
		The 2nd iteration	6.8909282035057489×10 -1
The 3rd iteration	2.4658347859604583×10 -1
		The 4th iteration	3.0359578642687003×10 -3

The 5th iteration	1.2278451094586332×10 -4
		The 6th iteration	5.2845269283793103×10 -5
The 7th iteration	6.0899512115445758×10 -7

Through 7 iteration, reach setting accuracy.

In the present invention, adopt differential correction algorithm pair $S_{\mathrm{BCM}}=\{{\mathrm{<figure-callout\; id="1"\; label="s"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">s</figure-callout>}}_1^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="2"\; label="s"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">s</figure-callout>}}_2^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="3"\; label="s"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">s</figure-callout>}}_3^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="4"\; label="s"\; filenames="HDA0000382880340000041.png,HDA0000382880340000081.png"\; state="\{\{state\}\}">s</figure-callout>}}_4^{\mathrm{BCM}}\};$ Revised, obtained the revised Lambert transfer orbit of first set ${\mathrm{DS}}_{\mathrm{BCM}}=\{{\mathrm{<figure-callout\; id="1"\; label="ds"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">ds</figure-callout>}}_1^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="2"\; label="ds"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">ds</figure-callout>}}_2^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="3"\; label="ds"\; filenames="HDA0000382880340000081.png"\; state="\{\{state\}\}">ds</figure-callout>}}_3^{\mathrm{BCM}},{\mathrm{<figure-callout\; id="4"\; label="ds"\; filenames="HDA0000382880340000041.png,HDA0000382880340000081.png"\; state="\{\{state\}\}">ds</figure-callout>}}_4^{\mathrm{BCM}}\},$ As shown in Fig. 5 A, Fig. 5 B: in figure, be heavy line after revising.The flight situation of spacecraft on the actual transfer track more approaches heavy line, more meets the characteristics of motion.Comparison diagram 5 and Fig. 5 A, wherein

for

correction,

for correction.Comparison diagram 5 and Fig. 5 B, wherein

for

correction,

for

correction.

(1) analyzing spacecraft is taken off to the intersection of Cycler track by ground

At first carry out the launch window analysis:

Red (earth equatorial plane) angle of cut is all monthly variation in vain, and the variation range during March 21 to 21 days March in 2006 in 2005 is at [0 °, η], and wherein the η value changes between 28.4 °～28.7 °.Xichang emission latitude is 27.9 °, and Wenchang emission latitude is 19 °, and both all exist launching opportunity the normal society face directly can be injected in space station, and the window of directly entering the orbit in Xichang is less than Wenchang; As shown in Figure 6, horizontal ordinate is the time, and ordinate is red (earth equatorial plane) angle of cut in vain.For Wenchang, the width of launch window 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, should select suitable launching time (1 day 1 time) so that right ascension of ascending node is captured in the normal society face.

The emission latitude in Xichang is with reference to " launching site, Xichang of " flue water outlet "---the sketch China four large satellite launching centres (1) " within 2007, delivered, the 17th page, author Xiao Bo.

The emission latitude of Wenchang is with reference to " the thunder and lightning environmental analysis of Wenchang, hainan rocketdrome " within 2012, delivered, the 183rd page, author Gao Yi etc.

Then carry out the intersection window analysis:

Resonance Cycler track is transferred to by parking orbit (being highly 100km) in space station, takes Huo Man branch mode or Lambert branch mode to carry out the analysis of Fuel Consumption and transfer time, as shown in Figure 6 A and 6 B.In Fig. 6 A, for fuel saving Huo Man shifts required speed increment, be 3192m/s, be 7 days transfer time.In Fig. 6 B, in order to take most fuel Huo Man, to shift required speed increment be 4093.5m/s, and be 66.5min transfer time.

Because reaching the apogean time, Fig. 6 A and Fig. 6 B (be about 7 days) about equally, and Fig. 6 A fuel is more saved, will be as nominal intersection track: note nominal intersection be the 1st day launching track, if because nominal intersection emission is incured loss through delay, select within the 2nd, 3,4,5,6,7 days, launched, required fuel increases to 4093.5m/s by 3192km/s successively; Note nominal intersection is the 7th day launching track, also can select, from day number scale reciprocal, within-6 ,-5 ,-4 ,-3 ,-2 ,-1 days, to be launched, and required fuel increases to 4093.5m/s by 3192km/s successively; Note nominal intersection is 1st～7 days launching tracks, and required fuel increases to 4093.5m/s by 3192km/s successively.As shown in Figure 6 C and 6 D shown in FIG., under ground moon barycenter inertial coordinates system O-XYZ, as shown in Figure 6 C, under ground moon barycenter rotating coordinate system O-xyz, launching track as shown in Figure 6 D for launching track for the launching track of continuous 7 days.

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

At first carry out the launch window analysis:

Red (moon equatorial plane) angle of cut is all monthly variation in vain, variation range during March 21 to 21 days March in 2006 in 2005 is at [0 °, ε], wherein the ε value changes between 6.64 °～6.86 °, as shown in Figure 7, horizontal ordinate is the time, and ordinate is red (moon equatorial plane) angle of cut in vain.Period of change is about 13.6 days.Landing point for latitude higher than 6.64 °, need to adjust orbit plane after lunar surface takes off, the resonance of being allowed for access Cycler track.Take the south poles landing point as example, and lunar surface takes off and enters 100km ring month SSO (Sun Synchronous Orbit); Resonance Cycler orbital acquisition process is as follows: the lifting apocynthion is to 20000km, and required speed increment is about 577m/s; Inclination correction is carried out in apocynthion, and required speed increment is about 287.3m/s.The lifting of apocynthion can be next step resonance Cycler orbital rendezvous and saves portion of energy.The moon revolution cycle is identical with the rotation period, only has every month twice chance right ascension of ascending node can be captured in the normal society face.

Then carry out the intersection window analysis:

Take intersection time and resonance Cycler orbital phase is variable, usings lunar surface parking orbit height 100km as constraint condition, and the search speed increment is less than or equal to the intersection window of 1500m/s.Geographic longitude depends on right ascension of ascending node, and according to " parking orbit height 100km " as the Search Results of constraint condition as shown in Figure 7 A, in Fig. 7 A, T1 point, T2 point, T3 point and T4 point mean to meet search condition and optimum intersection situation.Geographic longitude as shown in Figure 7 B can be distributed in [170 ° ,-135 °] ∪ [43 ° ,-5 °] ∪ [45 °, 108 °] ∪ [150 °, 165 °].Landing point in above-mentioned interval, be expected to directly enter in the normal society face after taking off by lunar surface; Landing point outside above-mentioned interval need to wait for that two weeks carries out right ascension of ascending node and catches.Therefore, if select the landing point in appointed area, after taking off, lunar surface enters in the normal society face: for the intersection of can taking off at any time in resonance Cycler orbital phase [0.05,0.25] ∪ [0.8,1].

In the present invention, adopt Gauss-global variable composite algorism to SM _lprocessed, obtained the second cover Lambert transfer orbit

to SM _ldisposal route with to SM _bCMdisposal route be identical.In order to distinguish explanation, add alphabetical H and distinguished on the relational expression of citation.

According to Lambert flight time theorem, the major semi-axis Hra of the transfer time on elliptical transfer orbit between any two points and elliptical transfer orbit, half past two footpath sum (HL ₁+ HL ₂) and central angle H θ relevant, be expressed as:

Ht _f=W(Hra,(HL ₁+HL ₂)，HR _P) (13)

If HL ₁with HL ₂sum is constant, and major semi-axis Hra is constant, origin endpoint E ₁with terminal E ₂between distance H R _pfor constant, from origin endpoint E ₁to terminal E ₂flight Ht transfer time _fit is also constant.Definite and selection poimt-to-point speed of elliptical transfer orbit is the key of Lambert problem, and it is described to following Gauss's problem: pursuit spacecraft P sets out and locates 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 ₁the initial velocity at some place is Hv ₁₀, spacecraft P is at elliptical transfer orbit E ₂the end speed at some place is Hv ₂₀.The target of this Gauss's problem is in order to solve initial position speed increment △ Hv ₁with terminal location speed increment Δ Hv ₂, and then the impulse force size applied of asking.

In the present invention, Gauss's problem can be tried to achieve by following transcendental equation group:

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

${\mathrm{Hv}}_{20}=\stackrel{\&CenterDot;}{k}\&times;{\mathrm{<figure-callout\; id="1"\; label="HL"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">HL</figure-callout>}}_1+\stackrel{\&CenterDot;}{g}\&times;{\mathrm{Hv}}_{10}---\left(15\right)$

Lagrange coefficient one $k=1-\frac{\mathrm{\&mu;}\&times;{\mathrm{<figure-callout\; id="2"\; label="HL"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">HL</figure-callout>}}_2}{{\mathrm{Hh}}^2}(1-\cos \mathrm{H\&theta;});$

Lagrange coefficient two $\stackrel{\&CenterDot;}{k}=\frac{\mathrm{\&mu;}}{\mathrm{Hh}}\frac{1-\cos \mathrm{H\&theta;}}{\sin \mathrm{H\&theta;}}[\frac{\mathrm{\&mu;}}{\mathrm{Hh}}(1-\cos \mathrm{H\&theta;})-\frac{1}{{\mathrm{HL}}_1}-\frac{1}{{\mathrm{HL}}_2}];$

Lagrange coefficient three $g=\frac{{\mathrm{<figure-callout\; id="1"\; label="HL"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">HL</figure-callout>}}_1\&times;{\mathrm{<figure-callout\; id="2"\; label="HL"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">HL</figure-callout>}}_2\&times;\sin \mathrm{H\&theta;}}{\mathrm{Hh}};$

Lagrange coefficient four

μ is time celestial body and primary body mass ratio, and in the BCM model, value is 0.000003003; H is the variable of spacecraft P on the elliptical transfer orbit at diverse location place.

In formula (14), formula (15), known HL ₁and Hv ₁₀, or HL ₂and Hv ₂₀just can determine the elliptical transfer orbit of spacecraft P operation.Obviously, once determined Lagrangian coefficient k, g, gauss's problem just can be readily solved.

In the present invention, for formula (14), formula (15), adopt the global variable algorithm to be solved.

Suppose that be 1 day the earth to homoclinic orbit transfer time, transfer orbit start point distance earth surface height 100km.Be 4 days homoclinic orbit to moon transfer time, and the transfer orbit terminal is apart from moonscape height 100km.According to disome Lambert problem solution, corresponding chummage Lambert track simulation result is as shown in Fig. 7 C.

Adopt differential correction algorithm pair

revised, obtained the revised Lambert transfer orbit of the second cover

the second cover is revised rear transfer orbit simulation result as shown in 7D.

According to HL ₁, HL ₂can calculate disome Lambert transfer orbit with H θ.Lambert becomes the transfer orbit of rail strategy to the out position intersection, and controlled quentity controlled variable is origin endpoint E ₁the impulse speed increment of position, its deflection is designated as H β; The differential correction algorithm will improve initial position speed increment △ Hv ₁reach Ht transfer time _f, to realize end point E ₂the intersection of position.Press close to the iterative initial value of true value, can guarantee the convergence of differential correction algorithm iterative process.For spacecraft P is directed into to target location L ₂(being theoretical terminal location vector), each iterative process all by orbit integration to HL ₂place, and Ht _fbe taken as this orbit integration time.Obviously, origin endpoint E ₁the orbital velocity Hv of position ₁₀changes delta Hv ₁to cause orbit integration time Ht _fvariation, be designated as Δ Ht _f.Examine or check the m time iteration (the front m-1 that once is designated as that iteration is m time, once be designated as m+1 after iteration m time, iteration be designated as for m time ought be last time), to end point E ₂the position vector of position is made single order Taylor and is launched, and can obtain:

$\mathrm{HL}({\mathrm{Ht}}_f+\mathrm{\&Delta;}{\mathrm{Ht}}_f,{\mathrm{<figure-callout\; id="10"\; label="Hv"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">Hv</figure-callout>}}_{10}+\mathrm{\&Delta;}{\mathrm{Hv}}_1)={\mathrm{<figure-callout\; id="2"\; label="HL"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">HL</figure-callout>}}_2+\frac{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="HL"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">HL</figure-callout>}}_2}{\&PartialD;{\mathrm{<figure-callout\; id="10"\; label="Hv"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">Hv</figure-callout>}}_{10}}\&times;\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="1"\; label="Hv"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">Hv</figure-callout>}}_1+{\mathrm{<figure-callout\; id="20"\; label="Hv"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">Hv</figure-callout>}}_{20}\&times;\mathrm{\&Delta;}{\mathrm{Ht}}_f---\left(16\right)$

HL (Ht _f+ Δ Ht _f, Hv ₁₀+ Δ Hv ₁) mean that spacecraft P terminal location actual vector is about transfer time, at elliptical transfer orbit E ₁function between the initial velocity at some place;

mean terminal location vector and elliptical transfer orbit E ₁partial derivative between the initial velocity at some place, be reduced to $\mathrm{HM}=\frac{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="HL"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">HL</figure-callout>}}_2}{\&PartialD;{\mathrm{<figure-callout\; id="10"\; label="Hv"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">Hv</figure-callout>}}_{10}};$

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{\&delta;}{\mathrm{<figure-callout\; id="2"\; label="HL"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">HL</figure-callout>}}_2={\mathrm{HL}}_2-{\stackrel{\&OverBar;}{\mathrm{HL}}}_2=\frac{\&PartialD;{\mathrm{<figure-callout\; id="2"\; label="HL"\; filenames="HDA0000382880340000021.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">HL</figure-callout>}}_2}{\&PartialD;{\mathrm{<figure-callout\; id="10"\; label="Hv"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">Hv</figure-callout>}}_{10}}\&times;\mathrm{\&Delta;}{\mathrm{Hv}}_1-{\mathrm{<figure-callout\; id="20"\; label="Hv"\; filenames="HDA0000382880340000011.png,HDA0000382880340000052.png"\; state="\{\{state\}\}">Hv</figure-callout>}}_{20}\&times;\mathrm{\&Delta;}{\mathrm{Ht}}_f---\left(18\right)$

the physical end position vector that means spacecraft P.

Might as well get

and L ₂there is identical horizontal ordinate component, at day ground rotating coordinate system O _s-X _sy _sz _sunder, and note terminal location vector and elliptical transfer orbit E ₁partial derivative between the initial velocity at some place

expansion (18) can obtain:

$\left[\begin{array}{c}0\\ \mathrm{\&delta;Hy}\\ \mathrm{\&delta;Hz}\end{array}\right]=-\left[\begin{array}{cccc}{\mathrm{HM}}_{11}& {\mathrm{<figure-callout\; id="12"\; label="HM"\; filenames="HDA0000382880340000052.png"\; state="\{\{state\}\}">HM</figure-callout>}}_{12}& {\mathrm{HM}}_{13}& H{v}_{20}^x\\ {\mathrm{HM}}_{21}& {\mathrm{<figure-callout\; id="22"\; label="HM"\; filenames="HDA0000382880340000052.png"\; state="\{\{state\}\}">HM</figure-callout>}}_{22}& {\mathrm{HM}}_{23}& H{v}_{20}^y\\ {\mathrm{HM}}_{31}& {\mathrm{HM}}_{32}& {\mathrm{HM}}_{33}& H{v}_{20}^z\end{array}\right]\&times;\left[\begin{array}{c}\mathrm{\&Delta;H}{v}_1^x\\ \mathrm{\&Delta;H}{v}_1^y\\ \mathrm{\&Delta;H}{v}_1^z\\ \mathrm{\&Delta;H}{t}_f\end{array}\right]---\left(19\right)$

δ Hy means a day ground rotating coordinate system O _s-X _sy _sz _sunder Y _sthe position difference of axle;

δ Hz means a day ground rotating coordinate system O _s-X _sy _sz _sunder Z _sthe position difference of axle;

HM ₁₁mean terminal location vector and elliptical transfer orbit E ₁the first row first row element of the partial derivative matrix between the initial velocity at some place;

HM ₁₂mean terminal location vector and elliptical transfer orbit E ₁the first row secondary series element of the partial derivative matrix between the initial velocity at some place;

HM ₁₃mean terminal location vector and elliptical transfer orbit E ₁the first row the 3rd column element of the partial derivative matrix between the initial velocity at some place;

HM ₂₁mean terminal location vector and elliptical transfer orbit E ₁the second row first row element of the partial derivative matrix between the initial velocity at some place;

HM ₂₂mean terminal location vector and elliptical transfer orbit E ₁the second row secondary series element of the partial derivative matrix between the initial velocity at some place;

HM ₂₃mean terminal location vector and elliptical transfer orbit E ₁the second row the 3rd column element of the partial derivative matrix between the initial velocity at some place;

HM ₃₁mean terminal location vector and elliptical transfer orbit E ₁the third line first row element of the partial derivative matrix between the initial velocity at some place;

HM ₃₂mean terminal location vector and elliptical transfer orbit E ₁the third line secondary series element of the partial derivative matrix between the initial velocity at some place;

HM ₃₃mean terminal location vector and elliptical transfer orbit E ₁the third line the 3rd column element of the partial derivative matrix between the initial velocity at some place;

mean terminal location E ₂end speed in place's is at X _sthe speed component of axle;

mean terminal location E ₂end speed in place's is at Y _sthe speed component of axle;

mean terminal location E ₂end speed in place's is at Z _sthe speed component of axle;

mean origin endpoint E ₁position initial velocity increment is at X _sthe speed component of axle;

mean origin endpoint E ₁position initial velocity increment is at Y _sthe speed component of axle;

mean origin endpoint E ₁position initial velocity increment is at Z _sthe speed component of axle;

△ Ht _fmean origin endpoint E ₁the orbital velocity v of position ₁₀increment △ v ₁the variation of the orbit integration time caused.

Order $\mathrm{HC}=\left[\begin{array}{cccc}{\mathrm{HM}}_{11}& {\mathrm{<figure-callout\; id="12"\; label="HM"\; filenames="HDA0000382880340000052.png"\; state="\{\{state\}\}">HM</figure-callout>}}_{12}& {\mathrm{HM}}_{13}& H{v}_{20}^x\\ {\mathrm{HM}}_{21}& {\mathrm{<figure-callout\; id="22"\; label="HM"\; filenames="HDA0000382880340000052.png"\; state="\{\{state\}\}">HM</figure-callout>}}_{22}& {\mathrm{HM}}_{23}& H{v}_{20}^y\\ {\mathrm{HM}}_{31}& {\mathrm{HM}}_{32}& {\mathrm{HM}}_{33}& H{v}_{20}^z\end{array}\right],$ Have:

$\left[\begin{array}{c}\mathrm{\&Delta;H}{v}_1^x\\ \mathrm{\&Delta;H}{v}_1^y\\ \mathrm{\&Delta;H}{v}_1^z\\ \mathrm{\&Delta;H}{t}_f\end{array}\right]=-{\left({\mathrm{HC}}^T\mathrm{HC}\right)}^{-1}{\mathrm{<figure-callout\; id="0"\; label="HC\; T"\; filenames="HDA0000382880340000012.png,HDA0000382880340000022.png"\; state="\{\{state\}\}">HC</figure-callout>}}^T\left[\begin{array}{c}0\\ \mathrm{\&delta;Hy}\\ \mathrm{\&delta;Hz}\end{array}\right]---\left(20\right)$

HC ^tthe inversion of representing matrix HC.

Therefore, with respect to the speed of a front m-1, last time the speed increment of m is characterized by:

$\mathrm{\&Delta;<figure-callout\; id="1"\; label="H\; v"\; filenames="HDA0000382880340000022.png,HDA0000382880340000041.png"\; state="\{\{state\}\}">H</figure-callout>}{v}_{}{}_1=\left[\begin{array}{ccc}\mathrm{\&Delta;H}{v}_1^x& \mathrm{\&Delta;H}{v}_1^y& \mathrm{\&Delta;H}{v}_1^z\end{array}\right]---\left(21\right)$

After setting accuracy, iterate, until meet precision.For moon parking orbit, take off to the transfer orbit of resonance type Cycler track, setting accuracy is 1 * 10 ^-15the iteration result is as shown in the table.

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

Through 10 iteration, reach setting accuracy.

Transfer orbit for earth parking orbit to chummage type Cycler track, setting accuracy is 1 * 10 ^-11the iteration result is as shown in the table.

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

Through 6 iteration, reach setting accuracy.

In the present invention, relatively based on resonance type Cycler track SM _bCMand chummage type Cycler track SM _lthe flight time of shifting task the ground moon (comprising the rail time that becomes), fuel consumption.Become the general speed momentum sign Fuel Consumption of rail with twice Lambert of each task.

1) resonance type Cycler track scheme:

By earth parking orbit to the total flight time of moon parking orbit, it is 6.014 days.

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

2) chummage type Cycler track scheme:

By earth parking orbit to the total flight time of moon parking orbit, it is 20.695 days.

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

By contrast total flight time, total fuel consumption, can find out, the ground moon transfer scheme based on the Cycler track, total flight time is longer, and total fuel consumption is less.This explanation, based on the Cycler track the ground moon transfer scheme imagination be feasible.Two class Cycler tracks have good periodicity.The cycle of resonance type Cycler track can be set, and spacecraft can, according to the task needs, in a planned way be travelled to and fro between ground between the moon, the transmission of carrying out goods and materials, personnel of rule.The cycle of chummage type Cycler track is also confirmable, and to have passed through the ground moon be LL ₁point, when completing ground month turnaround mission, the moon is LL over the ground ₁the scientific equipment that the some place arranges is changed and is overhauled.The transfer scheme flight time ground moon based on this two class Cycler track is longer, and total fuel consumption is less.Rationally utilize its periodically good advantage, can design ground month and come and go " 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 (8)

1. One kind the ground moon based on the Cycler track round task simulation system, it is characterized in that: this analogue system includes and builds resonance Cycler model trajectory (10), builds chummage Cycler model trajectory (20), multiple shooting method track correcting module (30) and Lambert transfer orbit acquisition module (40);

   Described structure resonance Cycler model trajectory (10) first aspect is according to CR3BP model construction M _cR3BPkinetic model; Second aspect adopts the BCM model to described M _cR3BPkinetic model is optimized, and obtains M _bCMkinetic model;

   Described structure chummage Cycler model trajectory (20) is according to BCM model construction M _lkinetic model;

   Described multiple shooting method track correcting module (30) first aspect adopts multiple shooting method respectively to M _bCMrevised, obtained revised resonance Cycler dynamics of orbits model DM _bCM; Second aspect is to DM _bCMadopt the fourth order Runge-Kutta method to obtain resonance Cycler track SM _bCM; The third aspect adopts multiple shooting method respectively to M _lrevised, obtained revised chummage Cycler dynamics of orbits model DM _l; Fourth aspect is to DM _ladopt the fourth order Runge-Kutta method to obtain chummage Cycler track SM _l;

   Described Lambert transfer orbit acquisition module (40) first aspect adopts Gauss-global variable composite algorism to SM _bCMprocessed, obtained first set Lambert transfer orbit $S_{\mathrm{BCM}}=\{s_1^{\mathrm{BCM}},s_2^{\mathrm{BCM}},s_3^{\mathrm{BCM}},s_4^{\mathrm{BCM}}\};$ Second aspect adopts differential correction algorithm pair

   

   revised, obtained the revised Lambert transfer orbit of first set ${\mathrm{DS}}_{\mathrm{BCM}}=\{{\mathrm{ds}}_1^{\mathrm{BCM}},{\mathrm{ds}}_2^{\mathrm{BCM}},{\mathrm{ds}}_3^{\mathrm{BCM}},{\mathrm{ds}}_4^{\mathrm{BCM}}\};$ The third aspect according to launch latitude and white red angle of cut Changing Pattern to SM _bCMcarry out the analysis of launch window and intersection window, the opportunity of obtaining Spacecraft Launch and entering the orbit; Fourth aspect adopts Gauss-global variable composite algorism to SM _lprocessed, obtained the second cover Lambert transfer orbit

   

   the 5th aspect adopts differential correction algorithm pair

   

   revised, obtained the revised Lambert transfer orbit of the second cover

   

   the 6th aspect is by comparatively Fuel Consumption, the total flight time of month two-way mission can be known SM _bCMand SM _ltrack advantage separately, for the ground moon two-way mission of the survey of deep space of low-thrust spacecraft provides the optimal design index.
2. The ground moon based on the Cycler track according to claim 1 round task simulation system, it is characterized in that: in described multiple shooting method track correcting module (30), at the M of one-period _bCMget the sampling point of 500 constant durations on model, any sampling point be designated as ini_b (i) (i=1 ..., 500); According to M _bCMmodel adopts fourth order Runge-Kutta method integration, the position and speed quantity of state that obtains any point be designated as in_a (i) (i=1 ..., 500); Ini_a (i :) is the matrix of 500 * 6; On each time interval, choose 100 sampling points, and sampling point is adopted to the ODE45 integration, the position and speed quantity of state φ of taking-up end point (ini_a (i :)); The position and speed quantity of state φ of note end point (ini_a (i :)) with the difference of the position and speed quantity of state ini_a of any point (i :), be F (i, :), i.e. F (i :)=φ (ini_a (i, :))-ini_a (i :); Application difference F (i :) makes the value of resolving and the actual computation value can be error free, in simulation process, should error be tapered to minimum as far as possible; Adopt Newton iteration method minimizing for above-mentioned difference F (i :).
3. The ground moon based on the Cycler track according to claim 1 round task simulation system, it is characterized in that: in described Lambert transfer orbit acquisition module (40), according to L ₁, L ₂can calculate disome Lambert transfer orbit with θ; Lambert becomes the transfer orbit of rail strategy to the out position intersection, and controlled quentity controlled variable is origin endpoint E ₁the impulse speed increment of position, its deflection is designated as β; The differential correction algorithm will improve initial position speed increment Δ v ₁reach t transfer time _f, to realize end point E ₂the intersection of position; For spacecraft P is directed into to target location L ₂, each iterative process all by orbit integration to L ₂place, and t _fbe taken as this orbit integration time; Obviously, origin endpoint E ₁the orbital velocity v of position ₁₀changes delta v ₁to cause orbit integration time t _fvariation, be designated as Δ t _f; Examine or check iteration the m time, to end point E ₂the position vector of position is made single order Taylor and is launched, and can obtain $L(t_f+\mathrm{\&Delta;}{t}_f,v_{10}+{\mathrm{\&Delta;v}}_1)=L_2+\frac{{\&PartialD;L}_2}{{\&PartialD;v}_{10}}\&times;{\mathrm{\&Delta;v}}_1+v_{20}\&times;{\mathrm{\&Delta;t}}_f;$ L(t _f+ Δ t _f, v ₁₀+ Δ v ₁) mean that spacecraft P terminal location actual vector is about transfer time, at elliptical transfer orbit E ₁function between the initial velocity at some place;

   

   mean terminal location vector and elliptical transfer orbit E ₁partial derivative between the initial velocity at some place, be reduced to

   

   Speed correction amount v ₁should make L (t _f+ Δ t _f, v ₁₀+ Δ v ₁)=L ₂set up, have ${\mathrm{\&delta;L}}_2=L_2-\stackrel{\&OverBar;}{L_2}=\frac{{\&PartialD;L}_2}{{\&PartialD;v}_{10}}\&times;\mathrm{\&Delta;}{v}_1-v_{20}\&times;\mathrm{\&Delta;}{t}_f;$

   

   the physical end position vector that means spacecraft P;

   Get

   

   and L ₂there is identical horizontal ordinate component, at day ground rotating coordinate system O _s-X _sy _sz _sunder, and note terminal location vector and elliptical transfer orbit E ₁partial derivative between the initial velocity at some place $M=\frac{{\&PartialD;L}_2}{{\&PartialD;v}_{10}},$ Have $\left[\begin{array}{c}0\\ \mathrm{\&delta;y}\\ \mathrm{\&delta;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]\&times;\left[\begin{array}{c}{\mathrm{\&Delta;v}}_1^x\\ {\mathrm{\&Delta;v}}_1^y\\ {\mathrm{\&Delta;v}}_1^z\\ {\mathrm{\&Delta;t}}_f\end{array}\right];$

   δ y means a day ground rotating coordinate system O _s-X _sy _sz _sunder Y _sthe position difference of axle;

   δ z means a day ground rotating coordinate system O _s-X _sy _sz _sunder Z _sthe position difference of axle;

   M ₁₁mean terminal location vector and elliptical transfer orbit E ₁the first row first row element of the partial derivative matrix between the initial velocity at some place;

   M ₁₂mean terminal location vector and elliptical transfer orbit E ₁the first row secondary series element of the partial derivative matrix between the initial velocity at some place;

   M ₁₃mean terminal location vector and elliptical transfer orbit E ₁the first row the 3rd column element of the partial derivative matrix between the initial velocity at some place;

   M ₂₁mean terminal location vector and elliptical transfer orbit E ₁the second row first row element of the partial derivative matrix between the initial velocity at some place;

   M ₂₂mean terminal location vector and elliptical transfer orbit E ₁the second row secondary series element of the partial derivative matrix between the initial velocity at some place;

   M ₂₃mean terminal location vector and elliptical transfer orbit E ₁the second row the 3rd column element of the partial derivative matrix between the initial velocity at some place;

   M ₃₁mean terminal location vector and elliptical transfer orbit E ₁the third line first row element of the partial derivative matrix between the initial velocity at some place;

   M ₃₂mean terminal location vector and elliptical transfer orbit E ₁the third line secondary series element of the partial derivative matrix between the initial velocity at some place;

   M ₃₃mean terminal location vector and elliptical transfer orbit E ₁the third line the 3rd column element of the partial derivative matrix between the initial velocity at some place;

   

   mean terminal location E ₂end speed in place's is at X _sthe speed component of axle;

   

   mean terminal location E ₂end speed in place's is at Y _sthe speed component of axle;

   

   mean terminal location E ₂end speed in place's is at Z _sthe speed component of axle;

   

   mean origin endpoint E ₁position initial velocity increment is at X _sthe speed component of axle;

   

   mean origin endpoint E ₁position initial velocity increment is at Y _sthe speed component of axle;

   

   mean origin endpoint E ₁position initial velocity increment is at Z _sthe speed component of axle;

   Δ t _fmean origin endpoint E ₁the orbital velocity v of position ₁₀increment Delta v ₁the variation of the orbit integration time caused;

   Order $C=\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],$ Have $\left[\begin{array}{c}{\mathrm{\&Delta;v}}_1^x\\ {\mathrm{\&Delta;v}}_1^y\\ {\mathrm{\&Delta;v}}_1^z\\ {\mathrm{\&Delta;t}}_f\end{array}\right]=-{\left(C^{T}C\right)}^{-1}{C}^T\left[\begin{array}{c}0\\ \mathrm{\&delta;y}\\ \mathrm{\&delta;z}\end{array}\right],$ C ^tthe inversion of representing matrix C; Therefore, with respect to the speed of a front m-1, last time the speed increment of m is characterized by ${\mathrm{\&Delta;v}}_1=\left[\begin{array}{ccc}{\mathrm{\&Delta;v}}_1^x& {\mathrm{\&Delta;v}}_1^y& {\mathrm{\&Delta;v}}_1^z\end{array}\right].$
4. The ground moon based on the Cycler track according to claim 1 round task simulation system, it is characterized in that: according to the CR3BP model Cycler dynamics of orbits model that obtains resonating $M_{\mathrm{CR}3\mathrm{BP}}=\left\{\begin{array}{c}\frac{\&PartialD;U}{\&PartialD;x}={\mathrm{HA}}^x-{2\mathrm{VA}}^y\\ \frac{\&PartialD;U}{\&PartialD;y}={\mathrm{HA}}^y+2{\mathrm{VA}}^x\\ \frac{\&PartialD;U}{\&PartialD;z}={\mathrm{HA}}^z\end{array}\right..$
5. The ground moon based on the Cycler track according to claim 1 round task simulation system, it is characterized in that: according to the BCM model, obtain chummage Cycler dynamics of orbits model $M_L=\left\{\begin{array}{c}{\mathrm{HC}}^{X_S}={2\mathrm{VC}}^{Y_S}+{\mathrm{xs}}_P-(1-{\mathrm{\&mu;}}_2)\frac{{\mathrm{xs}}_P+{\mathrm{\&mu;}}_2}{r_{\mathrm{PS}}^2}-\\ {\mathrm{\&mu;}}_2\frac{{\mathrm{xs}}_P-1+{\mathrm{\&mu;}}_2}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{xs}}_P-x_{P_2}}{r_{\mathrm{PM}}^3}\\ {\mathrm{HC}}^{Y_S}=-{2\mathrm{VC}}^{X_S}+{\mathrm{ys}}_P-(1-{\mathrm{\&mu;}}_2)\frac{{\mathrm{ys}}_P}{r_{\mathrm{PS}}^3}-\\ {\mathrm{\&mu;}}_2\frac{{\mathrm{ys}}_P}{r_{\mathrm{PE}}^3}-m_M\frac{{\mathrm{ys}}_P-y_{P_2}}{r_{\mathrm{PM}}^3}\\ {\mathrm{HC}}^{Z_S}=-(1-{\mathrm{\&mu;}}_2)\frac{{\mathrm{zs}}_P}{r_{\mathrm{PS}}^3}-{\mathrm{\&mu;}}_2\frac{{\mathrm{zs}}_P}{r_{\mathrm{PE}}^3}-m\frac{{\mathrm{zs}}_P}{r_{\mathrm{PM}}^3}\end{array}\right..$
6. The ground moon based on the Cycler track according to claim 1 round task simulation system, it is characterized in that: utilize launch latitude and white red angle of cut Changing Pattern gained, resonance Cycler track is transferred to by parking orbit in space station, take Huo Man branch mode or Lambert branch mode to carry out the analysis of Fuel Consumption and transfer time, it is 3192m/s that fuel saving Huo Man shifts required speed increment, and be 7 days transfer time; Taking most fuel Huo Man, to shift required speed increment be 4093.5m/s, and be 66.5min transfer time; Fuel saving Huo Man is shifted as nominal intersection track, if, because nominal intersection emission is incured loss through delay, select within the 2nd, 3,4,5,6,7 days, launched, required fuel increases to 4093.5m/s by 3192km/s successively; Note nominal intersection is the 7th day launching track, also can select, from day number scale reciprocal, within-6 ,-5 ,-4 ,-3 ,-2 ,-1 days, to be launched, and required fuel increases to 4093.5m/s by 3192km/s successively.
7. The ground moon based on the Cycler track according to claim 1 round task simulation system, it is characterized in that: utilize launch latitude and white red angle of cut Changing Pattern gained, take intersection time and resonance Cycler orbital phase is variable, using lunar surface parking orbit height 100km as constraint condition, and the search speed increment is less than or equal to the intersection window of 1500m/s; Geographic longitude depends on right ascension of ascending node, according to parking orbit height l00km, as constraint condition, can search for the optimum intersection situation that obtains; Landing point in above-mentioned interval, be expected to directly enter in the normal society face after taking off by lunar surface; Landing point outside above-mentioned interval need to wait for that two weeks carries out right ascension of ascending node and catches.
8. The ground moon based on the Cycler track according to claim 1 round task simulation system, it is characterized in that: relatively based on resonance type Cycler track SM _bCMand chummage type Cycler track SM _lthe flight time, the fuel consumption that shift task the ground moon; Become the general speed momentum sign Fuel Consumption of rail with twice Lambert of each task;

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

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

Images