CN106547991A - Along the disturbance gravitation reconstruction model optimization method of glide trajectories - Google Patents

Abstract

The invention discloses a kind of disturbance gravitation reconstruction model optimization method along glide trajectories, comprises the following steps：First, calculate change poles coordinate system trajectory three-dimensional envelope；Secondly, carry out change poles coordinate system universe spatial domain subdivision；Again, set up General Coordinate System disturbance gravitation universe reconstruction model；Then, set up the disturbance gravitation partial reconfiguration model along trajectory；Based on this, disturbance gravitation Fast approximating in flight course is set up；Afterwards, carry out EXPERIMENTAL DESIGN and build the agent model of reconstruction model；Finally, optimization method is set up for agent model.Carry out optimizing along the disturbance gravitation reconstruction model of glide trajectories according to the inventive method, the precision highest model under amount of storage least model and given amount of storage under given accuracy is required are required can be obtained, the requirement calculated on bullet in real time to memory data output, calculating speed and reconstruction accuracy is disclosure satisfy that, with preferable future in engineering applications.

Description

Along the disturbance gravitation reconstruction model optimization method of glide trajectories

Technical field

The invention belongs to vehicle dynamics models field, the more particularly to disturbance gravitation reconstruction model along glide trajectories is excellent Change method.

Background technology

Glide trajectories according to the present invention refer in particular to hypersonic glide vehicle comprising initial descending branch trajectory and gliding section In interior trajectory, its starting point is reentry point to trajectory, and terminal is that gliding Duan Zhongmo hands over to the next shift a little.Glide trajectories long-time is in closes on sky Between, affected significantly by earth disturbance gravitation.Calculating shows, the glide trajectories terminal location deviation caused by disturbance gravitation is up to several Ten kilometers, it is therefore necessary to the impact of disturbance gravitation is considered in ballistic computation and guidance are resolved.Existing disturbance gravitation calculating side Method generally has the characteristics of calculation scale is big, memory data output is big, the calculating time is long, it is impossible to meets Change the requirement with computational efficiency.To seek optimum reconstruction model parameter, a kind of disturbance gravitation reconstruct mould along glide trajectories is proposed Type optimization method.Method calculates change poles coordinate system trajectory three-dimensional envelope first；Next carries out change poles coordinate system universe spatial domain subdivision； General Coordinate System disturbance gravitation universe reconstruction model is set up again；Then set up the disturbance gravitation partial reconfiguration model along trajectory； Disturbance gravitation Fast approximating in flight course is set up based on this；Carry out EXPERIMENTAL DESIGN afterwards and build the agency of reconstruction model Model；Finally optimization method is set up for agent model.Method can obtain given accuracy require under amount of storage least model with And the precision highest model under given amount of storage requirement, with preferable future in engineering applications.

The content of the invention

The present invention is proposed a kind of based on agent model for the disturbance gravitation reconstruction model optimization problem along glide trajectories The disturbance gravitation reconstruction model optimization method along glide trajectories of reconstruction model, can obtain given accuracy require under amount of storage most Precision highest model under mini Mod and the requirement of given amount of storage, and with high precision, the feature that speed is fast, amount of storage is little.

The present invention is achieved through the following technical solutions above-mentioned purpose：

A kind of disturbance gravitation reconstruction model optimization method along glide trajectories, comprises the following steps：

The first step, sets up the mathematical model of trajectory envelope lateral width and vertical journey；

Second step, sets up change poles coordinate system；

3rd step, sets up vehicle dynamics model in change poles coordinate system；

4th step, calculates change poles coordinate system trajectory three-dimensional envelope；

5th step, change poles coordinate system universe spatial domain subdivision；

6th step, sets up General Coordinate System disturbance gravitation universe reconstruction model；

7th step, sets up the disturbance gravitation partial reconfiguration model along trajectory；

8th step, sets up；

9th step, obtains agent model training sample design conditions；

Tenth step, sets up the agent model of reconstruction model；

11st step, sets up the optimized algorithm of reconstruction model.

Specifically, the glide trajectories refer in particular to a part for hypersonic glide vehicle overall trajectory, and point of curve is Reentry point, ballistic impact are handed over to the next shift a little for middle end.

Above-mentioned steps are described further：

The first step, sets up the mathematical model of trajectory envelope lateral width and vertical journey

Trajectory reentry point is made to be I, its longitude is λ_I, reduced latitude be φ_I；Ballistic impact is T, and its longitude is λ_T, the earth's core latitude Spend for φ_T；Orthodrome plane is reentered as the plane of symmetry with what is determined by reentry point and drop point, is claimed along the directive direction plane of symmetry with left bullet Road is left side Maneuver Ballistic Trajectory, along the directive direction plane of symmetry with right trajectory as right side Maneuver Ballistic Trajectory；The maximum left side Maneuver Ballistic Trajectory of note away from The ultimate range of the plane of symmetry is left margin B_l, maximum ultimate range of the right side Maneuver Ballistic Trajectory away from the plane of symmetry of note is right margin B_r；

According to given aircraft reentry point flight status parameter, endgame flight status parameter, flight course constraint Condition and end conswtraint condition, are (λ for reentry point_I, φ_I0 °, 0 ° of)=(), drop point is (λ_ki, 0 °) (i=1,2,3,4 ...) Maneuver Ballistic Trajectory carries out ballistic computation on a large scale in low latitude；Remember that the vertical journey of the corresponding trajectory of each drop point is L_i, L is calculated by formula (1)_i,

L_i=R_earccos(sinφ_Isinφ_T+cosφ_Icosφ_Tcos(λ_T-λ_I)) (1)

Wherein, R_eFor earth radius；

Maximum left side Maneuver Ballistic Trajectory and maximum right side are obtained to Maneuver Ballistic Trajectory by changing reentry point initial velocity azimuth； Remember that the lateral envelope width of the corresponding trajectory of each drop point is W_io；W is calculated by formula (2)_io,

W_io=B_l-B_r (2)

Consider modeling error, take C_wTimes W_iO is actual lateral envelope width W_i, it is rightAnd W_i(i=1, 2,3,4 ...) it is fitted, the trajectory envelope lateral width and the mathematical model of vertical journey as shown in formula (3) can be set up, mould is determined Type coefficient a_i,

Second step, sets up change poles coordinate system

Introduce a change poles coordinate system；It is convenient for statement, useEach physical quantity in change poles coordinate system is represented, one is represented with X As each physical quantity in coordinate system；Change poles coordinate system is set up as follows：

1. define one orthodrome plane is reentered as change poles equatorial plane：1) situation about determining to impact point, will reenter What point and impact point the earth's core radius vector were constituted reenters orthodrome plane as change poles equatorial plane；2) for the undetermined feelings of impact point Condition, reenters orthodrome plane as the change poles equatoriat plane according to what reentry point position and Velocity Azimuth angle determined；

2. change poles coordinate system is defined based on change poles equatorial planeO_EFor the earth's core,Axle is sweared along reentry point the earth's core Footpath direction,Axle in the change poles equatoriat plane perpendicular toAxle points to impact point direction,Axle withAxle,Axle constitutes right-handed system；

3rd step, sets up vehicle dynamics model in change poles coordinate system

The glide vehicle kinetics equation with the time as independent variable is set up in change poles coordinate system, its state of flight amount is Longitude after change polesReduced latitudeFlight path yaw angleSpeedSpeed inclination angleWith the earth's core away from

Wherein, C_σ、C_θFor Corioli's acceleration item,WithFor aceleration of transportation item,

Wherein,

Wherein, ω_eFor earth rotation acceleration, λ_pAnd φ_pFor the longitude and reduced latitude of limit P after change poles, A_PFor The azimuth of P；

Defined according to change poles coordinate system, in General Coordinate System and change poles coordinate system, the earth's core is away from, local speed inclination angle and speed Definition it is consistent,

Definition

Wherein, ψ_fFor the azimuth of point F；

Determined in change poles coordinate system by λ in General Coordinate System and φWithExpression formula be,

By in change poles coordinate systemWithThe expression formula for determining λ and φ in General Coordinate System is,

Determined in change poles coordinate system by σ in General Coordinate SystemExpression formula be

Wherein,

4th step, calculates change poles coordinate system trajectory three-dimensional envelope

The coordinate transformation relation according to the 3rd step, can be according to trajectory reentry point longitude λ_I, reentry point reduced latitude φ_I, fall Point longitude λ_T, drop point reduced latitude φ_TIt is calculated change poles trajectory parameter；By calculatingNote Change poles ballistic impact longitude is

The vertical journey of change poles trajectory is calculated by formula (15)

Wherein, R_eFor earth radius；

WillSubstitution formula (16), can be calculated the lateral envelope width of change poles trajectory

Wherein, model coefficient a_iDetermined by the first step；

Change poles trajectory lateral envelope is described as length isWidth isRectangle, its border by formula (17) determine,

Wherein,For the lateral envelope east orientation lower boundary of change poles system,For the lateral envelope east orientation coboundary of change poles system；For The lateral envelope north orientation lower boundary of change poles system,For the lateral envelope north orientation coboundary of change poles system；

The lateral corresponding longitude range Δ λ of envelope and reduced latitude range delta φ determined by formula (18),

Consider motor-driven and deviation, with the C of ballistic ordinate scope_rAgain as the vertical scope of trajectory envelope

5th step, change poles coordinate system universe spatial domain subdivision

Note disturbance gravitation universe reconstruction model is Field models, and the spatial domain which determines is Ω；To ensure that first and last is prolonged Open up a little meaningful, change poles system Ω borders are calculated by formula (19),

Wherein,It is change poles system Ω days to lower boundary,It is change poles system Ω days to coboundary；For change poles system Ω east orientations Lower boundary,For change poles system Ω east orientations coboundary；For change poles system Ω north orientation lower boundaries,For change poles system Ω north orientations coboundary； Dr, d λ and d φ are respectively spatial domain subdivision interval；

The change poles system universe spatial domain Ω determined by formula (19) is uniformly cutd open to the interval of dr, east orientation d λ, north orientation d φ according to day It is divided into the subdomain Ω of q non-overlapping copies_e(e=1,2 ..., q), note mesh point coordinate is

6th step, sets up General Coordinate System disturbance gravitation universe reconstruction model

According to change poles system mesh point coordinateBy described in the 3rd step, coordinate transformation relation calculates general coordinate It is mesh point coordinate N (r_G,λ_G,φ_G)；

General Coordinate System grid node disturbance gravitation position T is calculated by spheric harmonic function method_G,

Wherein,

With T_GCorresponding gravitation as disturbs gravitational acceleration δ g, i.e.,

δ g=gradT_G (23)

Gravitational acceleration δ g are disturbed then in day northeast coordinate system O_EThree component δ g in-REN_R、δg_E、δg_NFor：

Storage General Coordinate System node location three-component and node disturbance gravitation three-component, complete to disturb the reconstruct of gravitation universe Model (Field models) builds；

7th step, sets up the disturbance gravitation partial reconfiguration model along trajectory

Remember that the disturbance gravitation partial reconfiguration model along trajectory is Channel models；

Trajectory planning is carried out according to aerial mission, reference trajectory is calculated；

Judge reference trajectory after three layers of continuation grid, the section after continuation grid is read from Field model datas Point position three-component and node disturbance gravitation three-component, complete to disturb gravitation partial reconfiguration model (Channel models) structure；

8th step, sets up

Make subdomain Ω_eAnd its spatial domain that adjacent mesh determines is continuation domain Ω '_e, haveNote Ω_eComprising node Number is r, Ω '_eComprising nodes be s, have s ＞ r；

Set up subdomain local curveilinear coordinates system, the cell node described in local coordinate system and calculating point coordinates；Subdomain Ω_e Local curveilinear coordinates by radius r_l=r_GThe sphere of+dr/2, longitude λ_l=λ_GThe meridian plane of+d λ/2, latitude φ_l=φ_G+d The intersection composition of the latitude circle of φ/2, origin l (r_l,λ_l,φ_l) be three intersections intersection point, origin local coordinate be l (0,0,0)；ξ、 Respectively along the day of origin l to, north orientation and east orientation, then in unit, local coordinate A (ξ, η, ζ) of arbitrfary point A (r, λ, φ) is for η, ζ,

Unit summit A_i(r_i,λ_i,φ_i) local coordinate A_i(ξ_i,η_i,ζ_i) be

Make continuation domain Ω '_eNode disturbance gravitation value is u_i；In subdomain local curveilinear coordinates system, node can be disturbed gravitation It is described as the function with regard to local coordinate,

u(ξ,η,ζ):R³→R (27)

In Ω_eOne approximate function U of upper construction u:Ω_e→ R, meets U (x_i)=u_i(i=1,2 ..., n)；In subdomain Ω_e On take trinary polynomial class

Its front t item is Interpolation-Radix-Function, and r ＜ t ＜ s, even

Order

Undetermined coefficient a is solved by formula (31),

A substitutions formula (29) for obtaining will be solved, you can gravitation value u is disturbed according to node_iSolve any point in continuation domain Disturbance gravitation value；

9th step, obtains agent model training sample design conditions

Based on total divisor design, orthogonal design or Latin hypercube design carry out experimental design, obtain it is different it is to be designed because Value of the element under varying level in scope of design；

Factor to be designed includes：

Reconstruction model stress and strain model interval dr, d λ, d φ described in (1) the 5th step；

Continuation domain node number s described in (2) the 8th steps；

Trinary polynomial class item number t described in (3) the 8th steps；

The scope of design of factor to be designed is setting value, and basis for selecting is reconstruct required precision；

Jing experimental designs, can obtain N_experiIndividual design factor combines Xⁱ=(dr, d λ, d φ, t, s), wherein i=1, 2,…,N_experi, N_experiFor setting value；

Tenth step, sets up the agent model of reconstruction model

According to the design factor combination X under the calculated varying level of the 9th stepⁱ=(dr, d λ, d φ, t, s) developments are examined Consider the glide trajectories emulation of disturbance gravitation；Ballistic Simulation of Underwater condition is determined by aerial mission；Disturbance gravitation is according to the first step to the 8th The described disturbance gravitation reconstructing method of step is calculated；

Jing Ballistic Simulation of Underwater, can obtain different designs factors combine XⁱThe corresponding Channel moulds of=(dr, d λ, d φ, t, s) Type nodes N_ch-nodeAnd glide trajectories terminal location deviation dD caused by disturbance gravitation；

With Xⁱ=(dr, d λ, d φ, t, s) is right as the input of training agent model, with Yⁱ=(N, dD) is used as training agency Output it is right, the functional relationship between acts on behalf of mould to output to set up description input by genetic algorithm or other approximating methods Type；

11st step, sets up the optimized algorithm of reconstruction model

Consider following two classes optimization problem：

(1) optimization problem one：Specified trajectory required precision, sets up amount of storage least model.

Object function：minN；

Variable to be designed：dr、dλ、dφ、t、s；

Constraints：dD∈[dD_min,dD_max]；

Variable-value scope to be designed：dr∈[dr_min,dr_max]、dλ∈[dλ_min,dλ_max]、dφ∈[dφ_min,d φ_max]、t∈[t_min,t_max]、s∈[s_min,s_max]；

(2) optimization problem two：Given amount of storage requirement, sets up reconstruction accuracy highest model.

Object function：mindD；

Variable to be designed：dr、dλ、dφ、t、s；

Constraints：N∈(0,N_max]；

Variable-value scope to be designed：dr∈[dr_min,dr_max]、dλ∈[dλ_min,dλ_max]、dφ∈[dφ_min,d φ_max]、t∈[t_min,t_max]、s∈[s_min,s_max]；

Choose a kind of global optimization method and solve the problems referred to above, the amount of storage that can finally set up under given accuracy is required is minimum Precision highest model under model and the requirement of given amount of storage.

So far, a kind of disturbance gravitation reconstruction model optimization side along glide trajectories can finally be set up through above-mentioned 11 step Method.Method can obtain the precision highest mould under the amount of storage least model and given amount of storage under given accuracy is required is required Type, and with high precision, the feature that speed is fast, amount of storage is little.Compared with existing Research foundation, method proposed by the present invention has Advantages below：

1) a kind of optimization method for gravitation reconstruction model is disturbed along glide trajectories is proposed first, and method can be rapidly completed Disturbance gravitation along any glide trajectories is calculated, with approximation accuracy height, the feature that calculating speed is fast, memory data output is little.

2) a kind of thinking for carrying out reconstruction model optimization design based on agent model is proposed, agent model is based on EXPERIMENTAL DESIGN Method and high accuracy approximating method build, and with higher fitting precision, the error of fitting of terminal location deviation can be controlled Within 0.5%, nodes are substantially error free.

3) optimization method set up can obtain the amount of storage least model under any given required precision, work as terminal location When deviation is defined to 5m, on the corresponding bullet of optimal models, memory data output is only 377, and can meet memory data output on bullet will Ask.

4) optimization method set up can obtain the precision highest model under any given amount of storage is required, when nodes are limited When fixed 400, the corresponding endgame position deviation of optimal models is only 4.030m, can meet reconstruction accuracy requirement.

Description of the drawings

Fig. 1 is the lateral envelope schematic diagram of different directive trajectories；

Fig. 2 is that east orientation is different indulges the lateral envelope schematic diagram of journey trajectory；

Fig. 3 is glide trajectories three-dimensional envelope schematic diagram；

Fig. 4 is that change poles coordinate system and state of flight amount define schematic diagram；

Fig. 5 is the position relationship schematic diagram of limit P after trajectory reentry point I and change poles；

Position relationship schematic diagrams of the Fig. 6 for change poles coordinate system flight path yaw angle；

Fig. 7 is lateral universe grid and local grid schematic diagram；

Fig. 8 is longitudinal universe grid and local grid schematic diagram；

Fig. 9 is three-dimensional subdomain and continuation domain schematic diagram；

Figure 10 is RBF neural network structure schematic diagram；

Specific embodiment

With reference to instantiation, the present invention is further illustrated：

Emulated based on CAV-H dummy vehicles, fundamental simulation condition setting is：1. quality m=908kg, the plane of reference Product Sm=0.48375m²；2. reentry point state parameter：Velocity magnitude V_I=6500m/, s speed inclination angle theta_I=0 °, height H_I= 80km；3. glide section termination state parameter：Velocity magnitude V_k=2500m/s, speed inclination angle theta_k=0 °, height H_k=30km, terminal Flight path driftage angular displacement is not more than ± 5 °；4. flight constraints condition：Maximum heat flow densityIt is maximum dynamic Pressure q_max=100kPa, maximum overload n_max=3g；5. flight journey is treated away from impact point at the end of terminal：S_togo=100km.

Simulation computer is configured to：Intel (R) Core (TM) i5-3470CPU 3.20GHz, inside save as 3.46GB.Software Environment is Window XP operating systems, and calculation procedure is based on VC++6.0 exploitations.

Which comprises the following steps that：

The first step, sets up the mathematical model of trajectory envelope lateral width and vertical journey

Trajectory reentry point is made to be I, its longitude is λ_I, reduced latitude be φ_I；Ballistic impact is T, and its longitude is λ_T, the earth's core latitude Spend for φ_T.Orthodrome plane is reentered as the plane of symmetry with what is determined by reentry point and drop point, is claimed along the directive direction plane of symmetry with left bullet Road is left side Maneuver Ballistic Trajectory, along the directive direction plane of symmetry with right trajectory as right side Maneuver Ballistic Trajectory；The maximum left side Maneuver Ballistic Trajectory of note away from The ultimate range of the plane of symmetry is left margin B_l, maximum ultimate range of the right side Maneuver Ballistic Trajectory away from the plane of symmetry of note is right margin B_r。

For the near space of different directives, Maneuver Ballistic Trajectory is emulated on a large scale, and its reentry point longitude is λ_I=0 °, again Access point reduced latitude is φ_I=0 °, drop point longitude λ_TWith drop point reduced latitude φ_TIt is shown in Table 1.

1 different directive ballistic impact coordinates of table

Directive	Positive north	Due east	Due south	Due west
					Longitude (°)	0	60	0	-60
Reduced latitude (°)	60	0	-60	0

It is calculated the lateral envelope schematic diagram of different directive trajectories and sees Fig. 1, the lateral envelope border of different directive trajectories is shown in Table 2.From Fig. 1 and Biao 2, due to the impact of earth rotation and aspherical factor, east orientation trajectory envelope is maximum, north-south trajectory bag Network takes second place, and west is minimum to trajectory envelope；North-south trajectory envelope has eastwards small range to offset.

The lateral envelope border of 2 different directive trajectories of table

For directive be due east, vertical Cheng Butong near space Maneuver Ballistic Trajectory is emulated on a large scale.Trajectory reentry point Jing Spend for λ_I=0 °, reentry point reduced latitude is φ_I=0；Drop point reduced latitude be 0 °, drop point longitude be respectively 50 °, 60 °, 70 ° and 80°.Remember that the vertical journey of the corresponding trajectory of each drop point is L_i, L is calculated by formula (32)_i,

L_i=R_earccos(sinφ_Isinφ_T+cosφ_Icosφ_Tcos(λ_T-λ_I)) (32)

Wherein, R_eFor earth radius.

It is calculated corresponding vertical journey L of different drop points_iIt is shown in Table 3.

The corresponding vertical journey of 3 different drop points of table

Drop point longitude (°)	50	60	70	80
					Vertical journey Li(km)	5560	6671	7784	8896

By the maximum left side Maneuver Ballistic Trajectory that changes reentry point initial velocity azimuth to obtain and determined by aircraft ability and Maximum right side is to Maneuver Ballistic Trajectory.The different lateral envelope schematic diagrams of journey trajectory of indulging are shown in Fig. 2.

Remember that the lateral envelope width of the corresponding trajectory of each drop point is W_io.W is calculated by formula (33)_io,

W_io=B_l-B_r (33)

Consider modeling error, take n=1.5, that is, take 1.5 times of W_ioFor actual lateral envelope width W_i.It is calculated different vertical The lateral envelope width of journey trajectory is shown in Table 4.

Table 4 is different to indulge the lateral envelope width of journey trajectory

Vertical journey Li(km)	Left margin Bl(km)	Right margin Br(km)	Lateral width Wio(km)	1.5 times of width Wi(km)
					5560	1538	1538	3076	4614
6671	1597	1597	3194	4791
					7784	1418	1340	2757	4136
8896	1028	797	1825	2738

To L_i(i=1,2,3,4 ...) and W_i(i=1,2,3,4 ...) is fitted, and can set up the bullet as shown in formula (34) The mathematical model of road envelope lateral width and vertical journey,

W=-1.064813L+0.000879L²-1.225047e^-7L³+4.607571e^-12L⁴ (34)

Maneuver Ballistic Trajectory three-dimensional envelope schematic diagram is shown in Fig. 3 near space on a large scale.

Second step, sets up change poles coordinate system

Introduce a change poles coordinate system.It is convenient for statement, useEach physical quantity in change poles coordinate system is represented, one is represented with X As each physical quantity in coordinate system.Change poles coordinate system is set up as follows：

1. define one orthodrome plane is reentered as change poles equatorial plane：1) situation about determining to impact point, will reenter What point and impact point the earth's core radius vector were constituted reenters orthodrome plane as change poles equatorial plane；2) for the undetermined feelings of impact point Condition, reenters orthodrome plane as the change poles equatoriat plane according to what reentry point position and Velocity Azimuth angle determined.

2. change poles coordinate system is defined based on change poles equatorial planeO_EFor the earth's core,Axle is sweared along reentry point the earth's core Footpath direction,Axle in the change poles equatoriat plane perpendicular toAxle points to impact point direction,Axle withAxle,Axle constitutes right-handed system. See Fig. 4.

3rd step, sets up vehicle dynamics model in change poles coordinate system

The glide vehicle kinetics equation with the time as independent variable is set up in change poles coordinate system, its state of flight amount is Longitude after change polesReduced latitudeFlight path yaw angleSpeedSpeed inclination angleWith the earth's core away from

Wherein, C_σ、C_θFor Corioli's acceleration item,WithFor aceleration of transportation item,

Wherein,

Wherein, ω_eFor earth rotation acceleration, λ_pAnd φ_pFor the longitude and reduced latitude of limit P after change poles, A_PFor The azimuth of P, sees Fig. 5.

Defined according to change poles coordinate system, in General Coordinate System and change poles coordinate system, the earth's core is away from, local speed inclination angle and speed Definition it is consistent,

Definition

Wherein, ψ_fFor the azimuth of point F.

Determined in change poles coordinate system by λ in General Coordinate System and φWithExpression formula be,

By in change poles coordinate systemWithThe expression formula for determining λ and φ in General Coordinate System is,

Determined in change poles coordinate system by σ in General Coordinate SystemExpression formula be that its location diagram is shown in Fig. 6,

Wherein,

4th step, calculates change poles coordinate system trajectory three-dimensional envelope

For reentry point longitude λ_I=0 °, reentry point reduced latitude φ_I=0 °, drop point longitude λ_T=30 °, drop point the earth's core latitude Degree φ_TManeuver Ballistic Trajectory is calculated=40 ° of near space on a large scale.The coordinate transformation relation according to the 3rd step, can calculate Obtain change poles coordinate system trajectory parameter

The vertical journey of change poles trajectory is calculated by formula (46)

Wherein, R_eFor earth radius.

By calculatingWillSubstitution formula (47), can be calculated the lateral envelope width of change poles trajectory

Wherein, model coefficient a_iDetermined by the first step.

By calculatingChange poles trajectory lateral envelope is described as length isWidth isSquare Shape, its border by formula (48) determine,

Wherein,For the lateral envelope east orientation lower boundary of change poles system,For the lateral envelope east orientation coboundary of change poles system；For The lateral envelope north orientation lower boundary of change poles system,For the lateral envelope north orientation coboundary of change poles system.

The lateral corresponding longitude range Δ λ of envelope and reduced latitude range delta φ determined by formula (49),

By can be calculated,Consider motor-driven and deviation, take

5th step, change poles coordinate system universe spatial domain subdivision

Note disturbance gravitation universe reconstruction model is Field models, and the spatial domain which determines is Ω.To ensure that first and last is prolonged Open up a little meaningful, change poles system Ω borders are calculated by formula (19),

Wherein,It is change poles system Ω days to lower boundary,It is change poles system Ω days to coboundary；For change poles system Ω east orientations Lower boundary,For change poles system Ω east orientations coboundary；For change poles system Ω north orientation lower boundaries,For change poles system Ω north orientations top Boundary.Dr, d λ and d φ are respectively spatial domain subdivision interval.

TakeCalculated change poles system universe spatial domain scope is shown in Table 5.

5 change poles system universe spatial domain scope of table

The change poles system universe spatial domain Ω determined by formula (50) is uniformly cutd open to the interval of dr, east orientation d λ, north orientation d φ according to day It is divided into the subdomain Ω of q non-overlapping copies_e(e=1,2 ..., q), note mesh point coordinate is

6th step, sets up General Coordinate System disturbance gravitation universe reconstruction model

According to change poles system mesh point coordinateThe coordinate transformation relation by described in the 3rd step calculates general seat Mark system mesh point coordinate N (r_G,λ_G,φ_G)。

General Coordinate System grid node disturbance gravitation position T is calculated by spheric harmonic function method_G,

Wherein,

With T_GCorresponding gravitation as disturbs gravitational acceleration δ g, i.e.,

δ g=gradT_G (54)

Gravitational acceleration δ g are disturbed then in day northeast coordinate system O_EThree component δ g in-REN_R、δg_E、δg_NFor：

Storage General Coordinate System node location three-component and node disturbance gravitation three-component, complete to disturb the reconstruct of gravitation universe Model (Field models) builds.

7th step, sets up the disturbance gravitation partial reconfiguration model along trajectory

Remember that the disturbance gravitation partial reconfiguration model along trajectory is Channel models.

Trajectory planning is carried out according to aerial mission, reference trajectory is calculated.

Judge reference trajectory after three layers of continuation grid, the section after continuation grid is read from Field model datas Point position three-component and node disturbance gravitation three-component, complete to disturb gravitation partial reconfiguration model (Channel models) structure.

For the Maneuver Ballistic Trajectory on a large scale of the near space for not considering no-fly zone described in the 4th step, its lateral universe grid with Local grid schematic diagram is shown in Fig. 7, and longitudinal universe grid is shown in Fig. 8 with local grid schematic diagram.From Fig. 7 and Fig. 8, local grid Among being included in global grid, and feasible trajectory can be covered, illustrate the universe reconstruction model and partial reconfiguration mould of this method foundation Type is feasible.

8th step, sets up

Make subdomain Ω_eAnd its spatial domain that adjacent mesh determines is continuation domain Ω '_e, haveTake Ω_eComprising node Number is r=8, Ω '_eComprising nodes be s=32.Three-dimensional subdomain and continuation domain schematic diagram are shown in Fig. 9.

Set up subdomain local curveilinear coordinates system, the cell node described in local coordinate system and calculating point coordinates.Subdomain Ω_e Local curveilinear coordinates by radius r_l=r_GThe sphere of+dr/2, longitude λ_l=λ_GThe meridian plane of+d λ/2, latitude φ_l=φ_G+d The intersection composition of the latitude circle of φ/2, origin l (r_l,λ_l,φ_l) be three intersections intersection point, origin local coordinate be l (0,0,0).ξ、 Respectively along the day of origin l to, north orientation and east orientation, then in unit, local coordinate A (ξ, η, ζ) of arbitrfary point A (r, λ, φ) is for η, ζ,

Unit summit A_i(r_i,λ_i,φ_i) local coordinate A_i(ξ_i,η_i,ζ_i) be

Make continuation domain Ω '_eNode disturbance gravitation value is u_i.In subdomain local curveilinear coordinates system, node can be disturbed gravitation It is described as the function with regard to local coordinate,

u(ξ,η,ζ):R³→R (58)

In Ω_eOne approximate function U of upper construction u:Ω_e→ R, meets U (x_i)=u_i(i=1,2 ..., n).In subdomain Ω_e On take trinary polynomial class front t items be Interpolation-Radix-Function, wherein t=20,

Then have,

The problems referred to above can be described as,

Wherein, undetermined coefficient a can be solved by formula (62),

Introduce Lagrange multiplier l=[l₁,l₂,…,l₈]^T, the Lagrangian shown in structural formula (63), then formula (63) minima minL is the solution of problem described by formula (62).

L (a, l)=(Ga-u)^T(Ga-u)+2(G_Ia-u_I)l (63)

According to the principle of optimality, minL can be determined by formula (64),

Write as matrix form, i.e.,

Wherein,

To improve the solution efficiency of formula (65), type function is introduced.A is noticed for square formation and reversible, order

28 rank Matrix for Inverse Problem then can be converted into the inversion problem of 20 ranks and 8 rank matrixes, amount of calculation is reduced.Solve Formula (65) can be obtained,

Main amount of calculation concentrates on matrix [D₁₁G^T D₁₂] solution on, and the value of the matrix is only and son that current point is located Domain and its continuation domain are relevant, therefore only just need to be updated after Fixed Initial Point exceeds subdomain.The a substitution formulas for obtaining will be solved (60), you can gravitation value u is disturbed according to node_iSolve the disturbance gravitation value of any point in continuation domain.

9th step, obtains agent model training sample design conditions

With t=20, illustrate as a example by the simplified situation of s=32, d λ=d φ.Variable to be designed be dr, d φ, its value Scope be 5km≤dr≤15km, 1 °≤d φ≤20 °.It is 50 to choose sample number, based on optimum Latin hypercube experimental design side Method carries out experimental design, and design result is as shown in table 6.

6 experimental design result of table

Tenth step, sets up the agent model of reconstruction model

According to the design factor combination X under the calculated varying level of the 9th stepⁱ=(dr, d φ) carries out consideration disturbance The glide trajectories emulation of gravitation.Ballistic Simulation of Underwater condition is：Vertical journey 6700km, 90 ° of initial orientation angle, reentry point and drop point are located at spy Mountainous region (λ₀=80 °, φ₀=42 °, λ_k=147 °, φ_k=19.3 °).Disturbance gravitation is according to the first step to the 8th step Disturbance gravitation reconstructing method is calculated.

Jing Ballistic Simulation of Underwater, can obtain different designs factors combine XⁱThe corresponding Channel model nodes of=(dr, d φ) Number N_ch-nodeAnd glide trajectories terminal location deviation dD caused by disturbance gravitation, it is shown in Table 7.

7 design factor of table combines corresponding nodes and endgame position deviation

Sequence number	dD(m)	N	Sequence number	dD(m)	N	Sequence number	dD(m)	N
									1	170.454	125	18	246.864	188	35	256.619	136
2	4.636	386	19	24.969	156	36	246.078	174
									3	39.530	180	20	70.219	160	37	233.817	125
4	133.851	181	21	131.519	137	38	141.783	196
									5	256.215	156	22	227.628	134	39	38.314	264
6	254.519	122	23	316.849	146	40	0.625	671
									7	55.071	139	24	25.626	221	41	38.089	159
8	514.103	170	25	7.107	200	42	12.182	169
									9	5.142	483	26	254.281	122	43	102.740	197
10	412.685	149	27	24.486	291	44	60.344	236
									11	498.275	138	28	141.783	154	45	24.842	196
12	0.428	540	29	5.260	159	46	477.121	113
									13	437.707	125	30	63.788	196	47	38.965	127
14	0.248	868	31	12.449	290	48	11.051	274
									15	387.344	172	32	292.103	182	49	78.377	146
16	416.843	135	33	380.286	147	50	7.484	299
									17	116.753	190	34	2.580	438	-	-	-

With Xⁱ=(dr, d λ, d φ, t, s) is right as the input of training agent model, with Yⁱ=(N, dD) is used as training agency Output it is right, by RBF neural set up description input to output between functional relationship agent model, see Figure 10.

4 groups of randomly selected design variables are input into, are compared by the result of calculation by agent model with archetype, Calculate the fitting precision of agent model.Agent model is shown in Table 8 to the error of fitting of endgame position deviation, and agent model is to section The error of fitting of points is shown in Table 9.As can be seen that agent model has higher fitting precision high, the fitting of terminal location deviation Error can be controlled within 0.5%, and nodes are substantially error free.

Error of fitting of 8 agent model of table to endgame position deviation

Error of fitting of 9 agent model of table to nodes

11st step, sets up the optimized algorithm of reconstruction model

Based on archipelago genetic algorithm for solving optimization problem, its major parameter is arranged and is shown in Table 10.

10 archipelago genetic algorithm parameter of table is arranged

Project	Numerical value
		Evolutionary generation	20
Island number	10
		Subgroup scale	10
Crossover probability	1.0
		Mutation probability	0.01
Mobility	0.01
		Migrate excessive algebraic quantity	5

The mathematical description of optimization problem one is：

Object function：minN；

Variable to be designed：dr、dφ；

Constraints：0m≤dD≤5m；

Variable-value scope to be designed：5km≤dr≤15km, 1 °≤d φ≤20 °.

Jing calculates the Bestgrid division for searching out meet the constraint condition for 875 times.The corresponding dr of optimum reconstruction model is 10.040km, d φ is 3.349 °, and terminal location deviation is 4.570m, and nodes are 377.Thus, establish given accuracy requirement Under amount of storage least model.

The mathematical description of optimization problem two is：

Object function：mindD；

Variable to be designed：dr、dφ；

Constraints：0 ＜ N≤400；

Variable-value scope to be designed：5km≤dr≤15km, 1 °≤d φ≤20 °.

Jing calculates the Bestgrid division for searching out meet the constraint condition for 846 times.The corresponding dr of optimum reconstruction model is 10.189km, d φ is 3.111 °, and terminal location deviation is 4.030m, and nodes are 399.Thus, establishing given amount of storage will Precision highest model under asking.

Summary simulation result can obtain to draw a conclusion：

1) reconstructing method can be rapidly completed the disturbance gravitation along glide trajectories and calculate, with computational accuracy is high and bullet on data The little feature of amount of storage.

2) agent model has higher fitting precision, and the error of fitting of terminal location deviation can be controlled within 0.5%, Nodes are substantially error free；

3) optimization method by setting up herein, it is possible to obtain amount of storage least model under any given required precision and Precision highest model under any given amount of storage requirement.When terminal location deviation is defined to 5m, the corresponding bullet of optimal models Upper memory data output is only 377, can meet the requirement of memory data output on bullet.When nodes limit 400, optimal models pair The endgame position deviation answered is only 4.030m, can meet reconstruction accuracy requirement.

4) method for proposing has that approximation accuracy is high, calculating speed is fast, the feature of wide accommodation, and possesses quick response The ability of aerial mission temporary changes.

The preferred embodiments of the present invention are the foregoing is only, the present invention is not limited to, for the skill of this area For art personnel, the present invention can have various modifications and variations.It is all within the spirit and principles in the present invention, made any repair Change, equivalent, improvement etc., should be included within the scope of the present invention.

Claims (3)

1. a kind of disturbance gravitation reconstruction model optimization method along glide trajectories, it is characterised in that：Comprise the following steps：

The first step, sets up the mathematical model of trajectory envelope lateral width and vertical journey；

Second step, sets up change poles coordinate system；

3rd step, sets up vehicle dynamics model in change poles coordinate system；

4th step, calculates change poles coordinate system trajectory three-dimensional envelope；

5th step, change poles coordinate system universe spatial domain subdivision；

6th step, sets up General Coordinate System disturbance gravitation universe reconstruction model；

7th step, sets up the disturbance gravitation partial reconfiguration model along trajectory；

8th step, sets up；

9th step, obtains agent model training sample design conditions；

Tenth step, sets up the agent model of reconstruction model；

11st step, sets up the optimized algorithm of reconstruction model.

2. the disturbance gravitation reconstruction model optimization method along glide trajectories according to claim 1, it is characterised in that：It is described Glide trajectories refer in particular to a part for hypersonic glide vehicle overall trajectory, and point of curve is reentry point, during ballistic impact is Hand over to the next shift a little at end.

3. the disturbance gravitation reconstruction model optimization method along glide trajectories according to claim 2, it is characterised in that：

The first step, sets up the mathematical model of trajectory envelope lateral width and vertical journey

Trajectory reentry point is made to be I, its longitude is λ_I, reduced latitude be φ_I；Ballistic impact is T, and its longitude is λ_T, reduced latitude be φ_T；Orthodrome plane is reentered as the plane of symmetry with what is determined by reentry point and drop point, along the directive direction plane of symmetry with left trajectory is called Left side Maneuver Ballistic Trajectory, along the directive direction plane of symmetry with right trajectory as right side Maneuver Ballistic Trajectory；The maximum left side Maneuver Ballistic Trajectory of note is away from symmetrical The ultimate range in face is left margin B_l, maximum ultimate range of the right side Maneuver Ballistic Trajectory away from the plane of symmetry of note is right margin B_r；

According to given aircraft reentry point flight status parameter, endgame flight status parameter, flight course constraints And end conswtraint condition, it is (λ for reentry point_I,φ_I0 °, 0 ° of)=(), drop point is (λ_ki, 0 °) (i=1,2,3,4 ...) it is low Empty Maneuver Ballistic Trajectory on a large scale carries out ballistic computation, remembers that the vertical journey of the corresponding trajectory of each drop point is L_i, L is calculated by formula (1)_i,

L_i=R_e arccos(sinφ_I sinφ_T+cosφ_I cosφ_T cos(λ_T-λ_I)) (1)

Wherein, R_eFor earth radius；

Maximum left side Maneuver Ballistic Trajectory and maximum right side are obtained to Maneuver Ballistic Trajectory by changing reentry point initial velocity azimuth, note is each The lateral envelope width of the corresponding trajectory of drop point is W_io, W is calculated by formula (2)_io,

W_io=B_l-B_r (2)

Consider modeling error, take C_wTimes W_ioFor actual lateral envelope width W_i, to L_i(i=1,2,3,4 ...) and W_i(i=1,2, 3,4 ...) it is fitted, the trajectory envelope lateral width and the mathematical model of vertical journey as shown in formula (3) can be set up, model is determined Coefficient a_i,

$W=\underset{i=0}{\overset{term-1}{\Σ}}{a}_i{L}^i---\left(3\right)$

Second step, sets up change poles coordinate system

A change poles coordinate system is introduced, is that statement is convenient, is usedEach physical quantity in change poles coordinate system is represented, general coordinate is represented with X Each physical quantity in system, sets up change poles coordinate system as follows：

1. define one orthodrome plane is reentered as change poles equatorial plane：1) situation determined by impact point, by reentry point and What impact point the earth's core radius vector was constituted reenters orthodrome plane as change poles equatorial plane；2) for the undetermined situation of impact point, Orthodrome plane is reentered as the change poles equatoriat plane according to what reentry point position and Velocity Azimuth angle determined；

2. change poles coordinate system is defined based on change poles equatorial planeO_EFor the earth's core,Axle is along reentry point the earth's core radius vector side To,Axle in the change poles equatoriat plane perpendicular toAxle points to impact point direction,Axle withAxle,Axle constitutes right-handed system；

3rd step, sets up vehicle dynamics model in change poles coordinate system

The glide vehicle kinetics equation with the time as independent variable is set up in change poles coordinate system, its state of flight amount is change poles Longitude afterwardsReduced latitudeFlight path yaw angleSpeedSpeed inclination angleWith the earth's core away from

$\left\{\begin{array}{c}\frac{d\hat{r}}{dt}=\hat{V}\sin \widehat{\mathrm{\θ}}\\ \frac{d\widehat{\mathrm{\λ}}}{dt}=\frac{\hat{V}\cos \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\σ}}}{\hat{r}\cos \widehat{\mathrm{\φ}}}\\ \frac{d\widehat{\mathrm{\φ}}}{dt}=\frac{\hat{V}\cos \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\σ}}}{\hat{r}}\\ \frac{d\widehat{\mathrm{\σ}}}{dt}=\frac{L\sin \mathrm{\υ}}{\hat{V}\cos \widehat{\mathrm{\θ}}}+\frac{\hat{V}}{\hat{r}}\cos \widehat{\mathrm{\θ}}\tan \widehat{\mathrm{\φ}}\sin \widehat{\mathrm{\σ}}-{\hat{g}}_{{\mathrm{\ω}}_e}\frac{\cos \widehat{\mathrm{\φ}}\sin \widehat{\mathrm{\σ}}}{\hat{V}\cos \widehat{\mathrm{\θ}}}+C_{\mathrm{\σ}}+{\tilde{C}}_{\mathrm{\σ}}\\ \frac{d\widehat{\mathrm{\θ}}}{dt}=\frac{L}{\hat{V}}\cos \mathrm{\υ}+\frac{\hat{V}}{\hat{r}}\cos \widehat{\mathrm{\θ}}+\frac{{\hat{g}}_r\cos \widehat{\mathrm{\θ}}}{\hat{V}}+\frac{{\hat{g}}_{{\mathrm{\ω}}_e}}{\hat{V}}(\cos \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\φ}}-\cos \widehat{\mathrm{\σ}}\sin \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\φ}})+C_{\mathrm{\θ}}+{\tilde{C}}_{\mathrm{\θ}}\\ \frac{d\hat{V}}{dt}=-D+{\hat{g}}_r\sin \widehat{\mathrm{\θ}}+{\hat{g}}_{{\mathrm{\ω}}_e}(\cos \widehat{\mathrm{\σ}}\cos \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\φ}}+\sin \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\φ}})+{\tilde{C}}_V\end{array}\right.---\left(4\right)$

Wherein, C_σ、C_θFor Corioli's acceleration item,WithFor aceleration of transportation item,

$\left\{\begin{array}{c}{C}_{\mathrm{\σ}}=(2{\mathrm{\ω}}_{ex}-2\tan \widehat{\mathrm{\θ}}({\mathrm{\ω}}_{ey}\sin \widehat{\mathrm{\σ}}+{\mathrm{\ω}}_{ez}\cos \widehat{\mathrm{\σ}}\left)\right)\\ {\tilde{C}}_{\mathrm{\σ}}=-\frac{\hat{r}}{\hat{V}\cos \widehat{\mathrm{\θ}}}({\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ey}\cos \widehat{\mathrm{\σ}}-{\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ez}\sin \widehat{\mathrm{\σ}})\\ C_{\mathrm{\θ}}=2({\mathrm{\ω}}_{ez}\sin \widehat{\mathrm{\σ}}-{\mathrm{\ω}}_{ey}\cos \widehat{\mathrm{\σ}})\\ {\tilde{C}}_{\mathrm{\θ}}=\frac{\hat{r}}{\hat{V}}\[{\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ey}\sin \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\σ}}+{\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ez}\sin \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\σ}}+({\mathrm{\ω}}_{ey}^2+{\mathrm{\ω}}_{ez}^2)\cos \widehat{\mathrm{\θ}}\]\\ {\tilde{C}}_V=\hat{r}\[-{\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ey}\cos \widehat{\mathrm{\θ}}\sin \widehat{\mathrm{\σ}}-{\mathrm{\ω}}_{ex}{\mathrm{\ω}}_{ez}\cos \widehat{\mathrm{\θ}}\cos \widehat{\mathrm{\σ}}+({\mathrm{\ω}}_{ey}^2+{\mathrm{\ω}}_{ez}^2)\sin \widehat{\mathrm{\θ}}\]\end{array}\right.---\left(5\right)$

Wherein,

$\left\{\begin{array}{c}{\mathrm{\ω}}_{ex}={\mathrm{\ω}}_e(\cos \widehat{\mathrm{\λ}}\cos \widehat{\mathrm{\φ}}{\mathrm{cos\φ}}_p\cos A_p+\sin \widehat{\mathrm{\λ}}\cos \widehat{\mathrm{\φ}}{\mathrm{cos\φ}}_p\sin A_p+\sin \widehat{\mathrm{\φ}}{\mathrm{sin\φ}}_p)\\ {\mathrm{\ω}}_{ey}={\mathrm{\ω}}_e(-\sin \widehat{\mathrm{\λ}}{\mathrm{cos\φ}}_p\cos A_p+\cos \widehat{\mathrm{\λ}}{\mathrm{cos\φ}}_p\sin A_p)\\ {\mathrm{\ω}}_{ez}={\mathrm{\ω}}_e(-\cos \widehat{\mathrm{\λ}}\sin \widehat{\mathrm{\φ}}{\mathrm{cos\φ}}_p\cos A_p-\sin \widehat{\mathrm{\λ}}\sin \widehat{\mathrm{\φ}}{\mathrm{cos\φ}}_p\sin A_p+\cos \widehat{\mathrm{\φ}}{\mathrm{sin\φ}}_p)\end{array}\right.---\left(6\right)$

Wherein, ω_eFor earth rotation acceleration, λ_pAnd φ_pFor the longitude and reduced latitude of limit P after change poles, A_PFor the side of P Parallactic angle；

Defined according to change poles coordinate system, General Coordinate System and the earth's core in change poles coordinate system are determined away from, locality speed inclination angle and speed It is adopted consistent,

$r=\hat{r},\mathrm{\θ}=\widehat{\mathrm{\θ}},V=\hat{V}---\left(7\right)$

Definition

$\begin{array}{c}\left[\begin{array}{ccc}{\mathrm{cos\φ}}_f{\mathrm{cos\λ}}_f& {\mathrm{cos\φ}}_f{\mathrm{sin\λ}}_f& {\mathrm{sin\φ}}_f\\ -{\mathrm{sin\ψ}}_f{\mathrm{sin\λ}}_f-{\mathrm{cos\ψ}}_f{\mathrm{sin\φ}}_f{\mathrm{cos\λ}}_f& {\mathrm{sin\ψ}}_f{\mathrm{cos\λ}}_f-{\mathrm{cos\ψ}}_f{\mathrm{sin\φ}}_f{\mathrm{sin\λ}}_f& {\mathrm{cos\ψ}}_f{\mathrm{cos\φ}}_f\\ {\mathrm{cos\ψ}}_f{\mathrm{sin\λ}}_f-{\mathrm{sin\ψ}}_f{\mathrm{sin\φ}}_f{\mathrm{cos\λ}}_f& -{\mathrm{cos\ψ}}_f{\mathrm{cos\λ}}_f-{\mathrm{sin\ψ}}_f{\mathrm{sin\φ}}_f{\mathrm{sin\λ}}_f& {\mathrm{sin\ψ}}_f{\mathrm{cos\φ}}_f\end{array}\right]\\ \stackrel{\mathrm{\Δ}}{=}\left[\begin{array}{ccc}{G}_{11}& G_{12}& G_{13}\\ G_{21}& G_{22}& G_{23}\\ G_{31}& G_{32}& G_{33}\end{array}\right]\end{array}---\left(8\right)$

Wherein, ψ_fFor the azimuth of point F；

$\left\{\begin{array}{c}{G}_{11}\cos \mathrm{\φ}\cos \mathrm{\λ}+G_{12}\cos \mathrm{\φ}\sin \mathrm{\λ}+G_{13}\sin \mathrm{\φ}\stackrel{\mathrm{\Δ}}{=}{k}_1\\ G_{21}\cos \mathrm{\φ}\cos \mathrm{\λ}+G_{22}\cos \mathrm{\φ}\sin \mathrm{\λ}+G_{23}\sin \mathrm{\φ}\stackrel{\mathrm{\Δ}}{=}{k}_2\\ G_{31}\cos \mathrm{\φ}\cos \mathrm{\λ}+G_{32}\cos \mathrm{\φ}\sin \mathrm{\λ}+G_{33}\sin \mathrm{\φ}\stackrel{\mathrm{\Δ}}{=}{k}_3\end{array}\right.---\left(9\right)$

$\left\{\begin{array}{c}{G}_{11}\cos \widehat{\mathrm{\φ}}\cos \widehat{\mathrm{\λ}}+G_{21}\cos \widehat{\mathrm{\φ}}\sin \widehat{\mathrm{\λ}}+G_{31}\sin \widehat{\mathrm{\φ}}\stackrel{\mathrm{\Δ}}{=}{\tilde{k}}_1\\ G_{12}\cos \widehat{\mathrm{\φ}}\cos \widehat{\mathrm{\λ}}+G_{22}\cos \widehat{\mathrm{\φ}}\sin \widehat{\mathrm{\λ}}+G_{32}\sin \widehat{\mathrm{\φ}}\stackrel{\mathrm{\Δ}}{=}{\tilde{k}}_2\\ G_{13}\cos \widehat{\mathrm{\φ}}\cos \widehat{\mathrm{\λ}}+G_{23}\cos \widehat{\mathrm{\φ}}\sin \widehat{\mathrm{\λ}}+G_{33}\sin \widehat{\mathrm{\φ}}\stackrel{\mathrm{\Δ}}{=}{\tilde{k}}_3\end{array}\right.---\left(10\right)$

Determined in change poles coordinate system by λ in General Coordinate System and φWithExpression formula be,

$\begin{array}{cc}\left\{\begin{array}{c}\cos \widehat{\mathrm{\λ}}=k_1/\sqrt{{k_1}^2+{k_2}^2}\\ \sin \widehat{\mathrm{\λ}}=k_2/\sqrt{{k_1}^2+{k_2}^2}\end{array}\right.& \left\{\begin{array}{c}\sin \widehat{\mathrm{\φ}}=k_3\\ \cos \widehat{\mathrm{\φ}}=\sqrt{{k_1}^2+{k_2}^2}\end{array}\right.\end{array}---\left(11\right)$

By in change poles coordinate systemWithThe expression formula for determining λ and φ in General Coordinate System is,

$\begin{array}{cc}\left\{\begin{array}{c}\cos \mathrm{\λ}={\tilde{k}}_1/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\\ \sin \mathrm{\λ}={\tilde{k}}_2/\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right.& \left\{\begin{array}{c}\sin \mathrm{\φ}={\tilde{k}}_3\\ \cos \mathrm{\φ}=\sqrt{{{\tilde{k}}_1}^2+{{\tilde{k}}_2}^2}\end{array}\right.\end{array}---\left(12\right)$

Determined in change poles coordinate system by σ in General Coordinate SystemExpression formula be

$\widehat{\mathrm{\σ}}=\mathrm{\σ}+\mathrm{\η}---\left(13\right)$

Wherein,

$\left\{\begin{array}{c}\sin \mathrm{\η}=\frac{\sin (\mathrm{\λ}-{\mathrm{\λ}}_P){\mathrm{cos\φ}}_P}{\cos \widehat{\mathrm{\φ}}}\\ \cos \mathrm{\η}=-\cos (A_P-\widehat{\mathrm{\λ}})\cos (\mathrm{\λ}-{\mathrm{\λ}}_P)+\sin (A_P-\widehat{\mathrm{\λ}})\sin (\mathrm{\λ}-{\mathrm{\λ}}_P){\mathrm{sin\φ}}_P\end{array}\right.---\left(14\right)$

4th step, calculates change poles coordinate system trajectory three-dimensional envelope

The coordinate transformation relation according to the 3rd step, can be according to trajectory reentry point longitude λ_I, reentry point reduced latitude φ_I, drop point Jing Degree λ_T, drop point reduced latitude φ_TIt is calculated change poles trajectory parameter；By calculatingNote change poles Ballistic impact longitude is

The vertical journey of change poles trajectory is calculated by formula (15)

$\hat{L}=R_e\arccos (sin{\widehat{\mathrm{\φ}}}_{I}sin{\widehat{\mathrm{\φ}}}_T+cos{\widehat{\mathrm{\φ}}}_{I}cos{\widehat{\mathrm{\φ}}}_{T}cos({\widehat{\mathrm{\λ}}}_T-{\widehat{\mathrm{\λ}}}_I\left)\right)---\left(15\right)$

Wherein, R_eFor earth radius；

WillSubstitution formula (16), can be calculated the lateral envelope width of change poles trajectory

$\hat{W}=\underset{i=0}{\overset{term-1}{\Σ}}{a}_i{\hat{L}}^i---\left(16\right)$

Wherein, model coefficient a_iDetermined by the first step；

Change poles trajectory lateral envelope is described as length isWidth isRectangle, its border by formula (17) determine,

$\left\{\begin{array}{c}{\hat{E}}_{\mathrm{\λ}1}={\mathrm{\λ}}_I\\ {\hat{E}}_{\mathrm{\λ}2}=\frac{180\hat{L}}{{\mathrm{\πR}}_e}\end{array}\right.,\left\{\begin{array}{c}{\hat{E}}_{\mathrm{\φ}1}=-\frac{90\hat{W}}{{\mathrm{\πR}}_e}\\ {\hat{E}}_{\mathrm{\φ}2}=\frac{90\hat{W}}{{\mathrm{\πR}}_e}\end{array}\right.---\left(17\right)$

Wherein,For the lateral envelope east orientation lower boundary of change poles system,For the lateral envelope east orientation coboundary of change poles system；For change poles It is lateral envelope north orientation lower boundary,For the lateral envelope north orientation coboundary of change poles system；

The lateral corresponding longitude range Δ λ of envelope and reduced latitude range delta φ determined by formula (18),

$\left\{\begin{array}{c}\mathrm{\Δ}\widehat{\mathrm{\λ}}=\frac{180\hat{L}}{{\mathrm{\πR}}_e}\\ \mathrm{\Δ}\widehat{\mathrm{\φ}}=\frac{180\hat{W}}{{\mathrm{\πR}}_e}\end{array}\right.---\left(18\right)$

Consider motor-driven and deviation, with the C of ballistic ordinate scope_rAgain as the vertical scope of trajectory envelope

5th step, change poles coordinate system universe spatial domain subdivision

Note disturbance gravitation universe reconstruction model is Field models, and the spatial domain which determines is Ω；To ensure first and last continuation point It is meaningful, change poles system Ω borders are calculated by formula (19),

$\left\{\begin{array}{c}{\hat{F}}_{r1}={\hat{r}}_I-\mathrm{\Δ}\hat{r}-2\·dr\\ {\hat{F}}_{r2}={\hat{r}}_I+2\·dr\end{array}\right.,\left\{\begin{array}{c}{\hat{F}}_{\mathrm{\λ}1}={\widehat{\mathrm{\λ}}}_I-2\·d\mathrm{\λ}\\ {\hat{F}}_{\mathrm{\λ}2}={\widehat{\mathrm{\λ}}}_I+\mathrm{\Δ}\widehat{\mathrm{\λ}}+2\·d\mathrm{\λ}\end{array}\right.,\left\{\begin{array}{c}{\hat{F}}_{\mathrm{\φ}1}={\widehat{\mathrm{\φ}}}_I-\mathrm{\Δ}\widehat{\mathrm{\φ}}/2\\ {\hat{F}}_{\mathrm{\φ}2}={\widehat{\mathrm{\φ}}}_I+\mathrm{\Δ}\widehat{\mathrm{\φ}}/2\end{array}\right.---\left(19\right)$

Wherein,It is change poles system Ω days to lower boundary,It is change poles system Ω days to coboundary；It is following for change poles system Ω east orientations Boundary,For change poles system Ω east orientations coboundary；For change poles system Ω north orientation lower boundaries,For change poles system Ω north orientations coboundary；dr、 D λ and d φ is respectively spatial domain subdivision interval；

According to day to the interval of dr, east orientation d λ, north orientation d φ by the uniform subdivisions of change poles system universe spatial domain Ω determined by formula (19) it is The subdomain Ω of q non-overlapping copies_e(e=1,2 ..., q), note mesh point coordinate is

$\mathrm{\Ω}=\underset{e=1}{\overset{q}{\∪}}{\mathrm{\Ω}}_e---\left(20\right)$

6th step, sets up General Coordinate System disturbance gravitation universe reconstruction model

According to change poles system mesh point coordinateBy described in the 3rd step, coordinate transformation relation calculates General Coordinate System net Lattice node coordinate N (r_G,λ_G,φ_G)；

General Coordinate System grid node disturbance gravitation position T is calculated by spheric harmonic function method_G,

$T_G=\frac{\mathrm{\μ}}{r_G}\underset{n=2}{\overset{s}{\Σ}}{\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}\cos {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\sin {\mathrm{m\λ}}_G){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)---\left(21\right)$

Wherein,

${\stackrel{\‾}{C}}_{nm}^{*}=\left\{\begin{array}{cc}0& \begin{array}{cc}n=2& m=0\end{array}\\ {\stackrel{\‾}{C}}_{nm}& else\end{array}\right.---\left(22\right)$

With T_GCorresponding gravitation as disturbs gravitational acceleration δ g, i.e.,

δ g=gradT_G (23)

Gravitational acceleration δ g are disturbed then in day northeast coordinate system O_EThree component δ g in-REN_R、δg_E、δg_NFor：

$\left\{\begin{array}{c}{\mathrm{\δg}}_R=\frac{\∂T}{\∂r_G}=-\frac{\mathrm{\μ}}{{r_G}^2}\underset{n=2}{\overset{s}{\Σ}}(n+1){\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}\cos {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\sin {\mathrm{m\λ}}_G){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)\\ {\mathrm{\δg}}_G=\frac{1}{r_G{\mathrm{cos\φ}}_G}\frac{\∂T}{\∂{\mathrm{\λ}}_G}=-\frac{1}{{\mathrm{cos\φ}}_G}\frac{\mathrm{\μ}}{{r_G}^2}\underset{n=2}{\overset{s}{\Σ}}{\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}m({\stackrel{\‾}{C}}_{nm}^{*}\sin {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\cos {\mathrm{m\λ}}_G){\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)\\ {\mathrm{\δg}}_G=\frac{1}{r_G}\frac{\∂T}{\∂{\mathrm{\φ}}_G}=\frac{\mathrm{\μ}}{{r_G}^2}\underset{n=2}{\overset{s}{\Σ}}{\left(\frac{a_e}{r_G}\right)}^n\underset{m=0}{\overset{n}{\Σ}}({\stackrel{\‾}{C}}_{nm}^{*}\cos {\mathrm{m\λ}}_G+{\stackrel{\‾}{S}}_{nm}\sin {\mathrm{m\λ}}_G)\frac{d{\stackrel{\‾}{P}}_{nm}\left({\mathrm{sin\φ}}_G\right)}{{\mathrm{d\φ}}_G}\end{array}\right.---\left(24\right)$

Storage General Coordinate System node location three-component and node disturbance gravitation three-component, complete to disturb gravitation universe reconstruction model (Field models) builds；

7th step, sets up the disturbance gravitation partial reconfiguration model along trajectory

Remember that the disturbance gravitation partial reconfiguration model along trajectory is Channel models；

Trajectory planning is carried out according to aerial mission, reference trajectory is calculated；

Judge reference trajectory after three layers of continuation grid, the node position after continuation grid is read from Field model datas Three-component and node disturbance gravitation three-component are put, completes to disturb gravitation partial reconfiguration model (Channel models) structure；

8th step, sets up

Make subdomain Ω_eAnd its spatial domain that adjacent mesh determines is continuation domain Ω '_e, haveNote Ω_eComprising nodes be R, Ω '_eComprising nodes be s, have s ＞ r；

Set up subdomain local curveilinear coordinates system, the cell node described in local coordinate system and calculating point coordinates；Subdomain Ω_eOffice Portion's curvilinear coordinate system is by radius r_l=r_GThe sphere of+dr/2, longitude λ_l=λ_GThe meridian plane of+d λ/2, latitude φ_l=φ_G+dφ/2 Latitude circle intersection composition, origin l (r_l,λ_l,φ_l) be three intersections intersection point, origin local coordinate be l (0,0,0)；ξ, η, ζ point Not along the day of origin l to, north orientation and east orientation, then in unit, local coordinate A (ξ, η, ζ) of arbitrfary point A (r, λ, φ) is,

$\left\{\begin{array}{c}\mathrm{\ξ}=r-r_l\\ \mathrm{\η}=rcos\mathrm{\φ}(\mathrm{\λ}-{\mathrm{\λ}}_l)\\ \mathrm{\ζ}=r(\mathrm{\φ}-{\mathrm{\φ}}_l)\end{array}\right.---\left(25\right)$

Unit summit A_i(r_i,λ_i,φ_i) local coordinate A_i(ξ_i,η_i,ζ_i) be

$\left\{\begin{array}{c}{\mathrm{\ξ}}_i=r_i-r_l\\ {\mathrm{\η}}_i=r_i{\mathrm{cos\φ}}_i({\mathrm{\λ}}_i-{\mathrm{\λ}}_l)\\ {\mathrm{\ζ}}_i=r_i({\mathrm{\φ}}_i-{\mathrm{\φ}}_l)\end{array}\right.---\left(26\right)$

Make continuation domain Ω '_eNode disturbance gravitation value is u_i；In subdomain local curveilinear coordinates system, can be by node disturbance gravitation description It is the function with regard to local coordinate,

u(ξ,η,ζ):R³→R (27)

In Ω_eOne approximate function U of upper construction u:Ω_e→ R, meets U (x_i)=u_i(i=1,2 ..., n)；In subdomain Ω_eOn take Trinary polynomial class

$\begin{array}{c}\left\{g_j(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})\right\}=\{1,\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ},{\mathrm{\ξ}}^2,{\mathrm{\η}}^2,{\mathrm{\ζ}}^2,\mathrm{\ξ}\mathrm{\η},\mathrm{\ξ}\mathrm{\ζ},\mathrm{\μ}\mathrm{\ζ},{\mathrm{\ξ}}^3,{\mathrm{\η}}^3,{\mathrm{\ζ}}^3,\\ {\mathrm{\ξ}}^2\mathrm{\η},{\mathrm{\ξ}}^2\mathrm{\ζ},{\mathrm{\ξ\η}}^2,{\mathrm{\ξ\ζ}}^2,{\mathrm{\η\ζ}}^2,{\mathrm{\η}}^2\mathrm{\ζ},\mathrm{\ξ}\mathrm{\η}\mathrm{\ζ},\mathrm{...}\}\end{array}---\left(28\right)$

Its front t item is Interpolation-Radix-Function, and r ＜ t ＜ s, even

$\begin{array}{cc}U(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})=\underset{j=1}{\overset{t}{\Σ}}{a}_j{g}_j(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})& (\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})\∈{\mathrm{\Ω}}_e\end{array}---\left(29\right)$

Order

$\left\{\begin{array}{ccc}G={\left\{g_j({\mathrm{\ξ}}_i,{\mathrm{\η}}_i,{\mathrm{\ζ}}_i)\right\}}_{ij}& i=r+1,r+2,...,s;& j=1,2,...t\\ G_I={\left\{g_j({\mathrm{\ξ}}_i,{\mathrm{\η}}_i,{\mathrm{\ζ}}_i)\right\}}_{ij}& i=1,2,...,r;& j=1,2,...t\\ u={\[u_{r+1},u_{r+2},...,u_s\]}^T& & \\ u_I={\[u_1,u_2,...,u_r\]}^T& & \\ a={\[a_1,a_2,...,a_t\]}^T& & \end{array}\right.---\left(30\right)$

Undetermined coefficient a is solved by formula (31),

$\left\{\begin{array}{c}minI\left(a\right)={(Ga-u)}^T(Ga-u)\\ \begin{array}{cc}s.t.& G_{I}a-u_I=0\end{array}\end{array}\right.---\left(31\right)$

A substitutions formula (29) for obtaining will be solved, you can gravitation value u is disturbed according to node_iSolve the disturbance of any point in continuation domain Gravitation value；

9th step, obtains agent model training sample design conditions

Experimental design is carried out based on total divisor design, orthogonal design or Latin hypercube design, different factors to be designed is obtained and is existed Value in scope of design under varying level；

Factor to be designed includes：

Reconstruction model stress and strain model interval dr, d λ, d φ described in (1) the 5th step；

Continuation domain node number s described in (2) the 8th steps；

Trinary polynomial class item number t described in (3) the 8th steps；

The scope of design of factor to be designed is setting value, and basis for selecting is reconstruct required precision；

Jing experimental designs, can obtain N_experiIndividual design factor combines Xⁱ=(dr, d λ, d φ, t, s), wherein i=1,2 ..., N_experi, N_experiFor setting value；

Tenth step, sets up the agent model of reconstruction model

According to the design factor combination X under the calculated varying level of the 9th stepⁱ=(dr, d λ, d φ, t, s) carries out consideration and disturbs The glide trajectories emulation of dynamic gravitation；Ballistic Simulation of Underwater condition is determined by aerial mission；Disturbance gravitation is according to the first step to the 8th step institute The disturbance gravitation reconstructing method stated is calculated；

Jing Ballistic Simulation of Underwater, can obtain different designs factors combine XⁱThe corresponding Channel models sections of=(dr, d λ, d φ, t, s) Points N_ch-nodeAnd glide trajectories terminal location deviation dD caused by disturbance gravitation；

With Xⁱ=(dr, d λ, d φ, t, s) is right as the input of training agent model, with Yⁱ=(N, dD) is defeated as training agency's It is right to go out, and sets up agent model of the description input to output functional relationship between by genetic algorithm or other approximating methods；

11st step, sets up the optimized algorithm of reconstruction model

Consider following two classes optimization problem：

(1) optimization problem one：Specified trajectory required precision, sets up amount of storage least model；

Object function：minN；

Variable to be designed：dr、dλ、dφ、t、s；

Constraints：dD∈[dD_min,dD_max]；

Variable-value scope to be designed：dr∈[dr_min,dr_max]、dλ∈[dλ_min,dλ_max]、dφ∈[dφ_min,dφ_max]、t∈ [t_min,t_max]、s∈[s_min,s_max]；

(2) optimization problem two：Given amount of storage requirement, sets up reconstruction accuracy highest model；

Object function：mindD；

Variable to be designed：dr、dλ、dφ、t、s；

Constraints：N∈(0,N_max]；

Variable-value scope to be designed：dr∈[dr_min,dr_max]、dλ∈[dλ_min,dλ_max]、dφ∈[dφ_min,dφ_max]、t∈ [t_min,t_max]、s∈[s_min,s_max]；

Choose a kind of global optimization method and solve above-mentioned optimization problem, the amount of storage that can finally set up under given accuracy is required is minimum Precision highest model under model and the requirement of given amount of storage.