Method for designing flight route of gliding section of high-speed gliding aircraft based on non-cooperative measurement data modeling

Abstract

The invention provides a high-speed gliding aircraft gliding section route design method based on non-cooperative measurement data modeling. The method comprises the steps that aircraft maneuvering characteristics are extracted from non-cooperative measurement data, and a model maneuvering characteristic database is established; establishing a dynamic model for the gliding section, solving an acceleration vector by using the maneuvering characteristic database and a state vector measured in real time, and performing integration to generate a flight route of the aircraft; designing a turning curvature sequence and a turning angle threshold sequence, obtaining a plurality of flight paths, and obtaining a reachable area of the aircraft; and generating a whole course route by taking the gliding section starting point as an initial value, taking the arrival point as an end point and passing through the passing points. According to the method, complex parameters such as the attack angle, the speed and the control surface angle of the aircraft and a real control equation do not need to be obtained, the acceleration of the aircraft in the gliding process can be solved under the condition that only part of the position and speed parameters of the aircraft are actually measured, and a starting point, an arrival point and a passing point are autonomously selected; and the controllable planning design of the whole flight route of the gliding section is realized.

Description

Technical Field

[0001] The invention belongs to the technical field of high-speed gliding aircraft route design, and in particular relates to a high-speed gliding aircraft gliding segment route design method based on non-cooperative measurement data modeling. Background Art

[0002] High-speed gliders glide at high Mach numbers in the atmosphere, with high speed and strong maneuverability, so they have strong penetration capabilities and are increasingly valued. Research on the glide route design of high-speed gliders is of great significance for the research on high-speed target attack and defense confrontation.

[0003] At present, the design of high-speed gliding aircraft routes is mainly based on the integral solution of control equations and dynamic equations. These methods require the mapping relationship between the aircraft state and aerodynamic parameters, the mapping relationship between the aircraft state and control torque, and the control force boundary conditions of the aircraft. It is also necessary to obtain the basic state parameters of the aircraft such as angle of attack, roll angle, yaw angle, speed, and control surface angle, and to build an accurate navigation control model. The high dependence on these complex parameters and models means that only the aircraft design department that designs the aircraft can realize the route design of high-speed gliding aircraft. Summary of the invention

[0004] The purpose of the present invention is to solve the problem in the prior art that the route design of a high-speed gliding aircraft is highly dependent on complex parameters and equations, resulting in that only the aircraft design department can implement the route design of the aircraft, and a method for designing the gliding segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling is provided. The method proposes a new high-speed gliding aircraft dynamics modeling method, and designs a route solution method from the gliding starting point to the waypoint, from the waypoint to the waypoint, and from the waypoint to the arrival point. Only the starting point, the arrival point, and the waypoint need to be selected, and the active route design of the gliding segment of the gliding aircraft can be realized without relying on the model design department, which is of great significance to the research on the attack and defense confrontation of high-speed gliding aircraft.

[0005] To achieve the above purpose, the technical solution provided by the present invention is:

[0006] A method for designing a glide segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling is characterized in that it includes the following steps:

[0007] Step 1, extracting high-speed gliding aircraft maneuvering characteristics from non-cooperative measurement data and establishing a model maneuvering characteristics database, includes the following sub-steps:

[0008] Step 1.1, extract maneuvering features from the single full-flight track data of the high-speed glider aircraft, the maneuvering features include the mass-to-drag ratio g M / D, lift-to-drag ratio g L / D , normal overload g L ,include:

[0009] Step 1.1.1, calculate the normal overload g L : Where g0 is the standard sea level gravity acceleration, a AL is the lift acceleration vector;

[0010] Step 1.1.2, calculate lift-to-drag ratio g L / D : In the formula, a AD is the drag acceleration;

[0011] Step 1.1.3, calculate the mass-to-resistance ratio g M / D : In the formula, ρ r is the atmospheric density at the location of the aircraft, v a is the aircraft airspeed vector;

[0012] Step 1.1.4, calculate the distance d between the aircraft and the target:

[0013]

[0014] In the formula, r F is the position vector of the spacecraft's Earth-centered fixed system, r F,Target is the target geocentric fixed position mass;

[0015] Step 1.2, establish a model maneuverability feature database, including:

[0016] Step 1.2.1, establish a characteristic parameter sequence, including:

[0017] Mass-to-resistance ratio maximum sequence C M / Dmax,Ma =[c M / Dmax,m1 ,…,c M / Dmax,mn ]、

[0018] Mass-to-resistance ratio minimum sequence C M / Dmin,Ma =[c M / Dmin,m1 ,…,c M / Dmin,mn ]、

[0019] Lift-to-drag ratio maximum sequence C L / Dmax,Ma =[c L / Dmax,m1 ,…,c L / Dmax,mn ]、

[0020] Lift-to-drag ratio minimum value sequence C L / Dmin,Ma =[c L / Dmin,m1 ,…,c / L / Dmin,mn ]、

[0021] Normal overload maximum sequence GLmax,Ma =[g Lmax,m1 ,…,g Lmax,mn ]、

[0022] Normal overload minimum sequence G Lmin,Ma =[g Lmin,m1 ,…,g Lmin,mn ]、

[0023] And the average flight altitude sequence H at the end of the flight mean,D =[h mean,d1 ,…,h mean,dn ],

[0024] Where Ma={m1,…,m n} is the Mach number sequence, D={d1,…,d n} is the remaining flight sequence;

[0025] Step 1.2.2, use the linear interpolation method to solve the maximum value c of the mass-to-drag ratio corresponding to the Mach number m M / Dmax,m , minimum value of mass-to-resistance ratio c M / Dmin,m , lift-to-drag ratio maximum value c L / Dmax,m , lift-to-drag ratio minimum value c L / Dmin,m , Normal overload maximum value g Lmax,m , minimum normal overload g Lmin,m , and the average flight altitude h corresponding to the remaining flight distance d mean,d ;

[0026] Step 2, establishing a dynamic model for the glide segment of the high-speed glider, using the established model maneuvering characteristic database and the real-time measured state vector to solve and integrate the acceleration vector to generate the flight route of the aircraft, including the following sub-steps:

[0027] Step 2.1: Solve for the angle θ between the aircraft speed and the horizontal direction a :

[0028]

[0029] In the formula, u r is the vertical upward direction vector of the aircraft's location;

[0030] Step 2.2, establishing a dynamic model and solving the acceleration vector, includes the following steps:

[0031] Step 2.2.1, aircraft motion state determination, including:

[0032] Step 2.2.1.1, calculate the angle δ between the current velocity vector and the velocity vector at the time of extrapolation initialization;

[0033] Step 2.2.1.2, calculate the centripetal acceleration L required for the aircraft to fly at a curvature C in the horizontal direction j ;

[0034] Step 2.2.1.3, determine the motion state of the aircraft, when θ I <θ ST When θ ST <θ I ≤0, the aircraft is judged to be in a balanced gliding state. I >0, the aircraft is judged to be in a ballistic flight state before entering a balanced gliding state, where θ ST is the angle threshold between the inertial system velocity and the horizontal direction of the aircraft in the balanced gliding state, θ I is the angle between the inertial system velocity and the horizontal direction of the aircraft at the current moment;

[0035] Step 2.2.2, calculate the lift acceleration L in the pull-up state k :

[0036] L k =a GKJ +g Lmax,m g0,

[0037] In the formula, a GKJ is the component of gravity acceleration on the plane where the lift vector of the aircraft is located, a GKJ =||a G ||cosθ a , a G is the acceleration due to gravity, g Lmax,m is the maximum overload of the aircraft at the current Mach number;

[0038] Step 2.2.3: Calculate the lift acceleration L in the equilibrium gliding state k ;

[0039] Step 2.2.4, calculate the lift acceleration L of ballistic flight k ;

[0040] Step 2.2.5, determine whether flight altitude control is required: calculate the remaining distance d between the aircraft and the target, and calculate the current flight altitude h t , solve the corresponding average flight height h mean,d , if the current flight altitude is h t <h mean,d , then no flight altitude control is performed and go to step 2.2.6; otherwise, the following flight altitude control is performed:

[0041] Calculate the lift acceleration L on the trajectory surface k , and update the elevation control integral correction value hint , the value of i:

[0042]

[0043] In the formula, c A 、c B 、c c h is the adjustment coefficient with a fixed value in advance. i-1 is the aircraft altitude at the last iteration, h i is the current altitude of the aircraft;

[0044] Step 2.2.6, synthesizing and correcting the lift acceleration to the total lift acceleration L, includes:

[0045] Step 2.2.6.1, calculate the maximum lift acceleration a that can be obtained at the current Mach number m of the aircraft in its current state ALmax ,m:

[0046] a AL max,m =c L / D max,m ·a AD max,m ,

[0047] In the formula, c M / D min,m is the minimum mass-to-drag ratio at the current Mach number, c L / D max,m is the maximum lift-to-drag ratio at the current Mach number;

[0048] Step 2.2.6.2, if L k >a AL max,m , then let: L k =a AL max,m , L j =0; if L k ≤a AL max,m , then calculate the maximum lift acceleration L that can be used for horizontal maneuvers j max :

[0049]

[0050] Step 2.2.6.3, determine the current lift acceleration L j , if |L j |>L jmax , then let L j =D j L jmax , D j is the lift direction sign, otherwise L j No change;

[0051] Step 2.2.6.4, calculate the resultant force L on the projectile perpendicular to the airspeed direction ki :

[0052]

[0053] Step 2.2.6.5, in L ki >g L max,m When L k , L j :

[0054]

[0055]

[0056] Step 2.2.6.6, calculate the total lift acceleration L:

[0057] Step 2.2.7, calculate the drag acceleration a AD :

[0058] Step 2.2.7.1, calculate the minimum drag acceleration a at the current Mach number AD min,m and lift acceleration a ADmin,m :

[0059] a ALmin,m =c L / Dmin,m ·a ADmin,m ,

[0060] In the formula, c M / Dmax,m is the maximum mass-to-drag ratio at the current Mach number, c L / Dmin,m is the minimum lift-to-drag ratio at the current Mach number;

[0061] Step 2.2.7.2, calculate the drag acceleration a AD :

[0062]

[0063] Step 2.2.8, acceleration vector synthesis and integration, includes:

[0064] Step 2.2.8.1, calculate the aircraft drag acceleration vector a AD and the lift acceleration vector a AL :

[0065] a AL =L k ·u ka +L j ·u ja ,

[0066] In the formula, u Ka is the direction vector perpendicular to the airspeed vector and coplanar with the gravity acceleration vector and the airspeed vector, u Jais the direction vector perpendicular to the airspeed vector and gravity acceleration vector;

[0067] Step 2.2.8.2, for Integrate and generate the aircraft trajectory, and terminate the iteration when the iteration termination condition is met: In the formula, v I is the velocity vector;

[0068] Step 3, designing the aircraft turning curvature sequence and the aircraft turning angle threshold sequence and bringing them into step 2 to obtain multiple flight tracks, thereby obtaining the aircraft reachable area, including the following sub-steps:

[0069] Step 3.1, design the curvature sequence C of the aircraft:

[0070] C=[c1,…,c n＝2m+3 ],m∈N + ,

[0071]

[0072] Where, d cr,j is the intermediate variable function, the subscript j is the independent variable of the function, c r2 is the maximum turning radius, c r1 is the minimum turning radius;

[0073] Step 3.2, design the aircraft turning angle threshold sequence δ threshold :

[0074] δ threshold =[δ1,…,δ n＝2m+3 ],m∈N + ,

[0075] Step 3.3: Combine the turning curvature and the turning angle threshold K i =[c i ,δ i ] T , i=1…n are successively brought into step 2, and the state vector s at a certain moment is I,t As the initial point of extrapolation, each flight route is integrated and calculated. When the iteration termination condition is met, the iteration is terminated until all routes are calculated, thereby generating n flight routes, which are connected in sequence with the endpoints of adjacent routes and the two outermost flight routes to form a closed envelope area.

[0076] The ground projection area of ​​the envelope area is the area accessible by the aircraft;

[0077] Step 4, taking the starting point of the gliding segment as the initial value and the arrival point as the end point, and passing through a number of passing points, generates the entire gliding segment route, including the following sub-steps:

[0078] Step 4.1, construct the starting point and the arrival point, including:

[0079] Step 4.1.1, select the starting point according to the measured data: s I,0 =[r I,0 ,v I,0 ] T , where s I,0 is the inertial system state vector at the start point of the aircraft glide segment, r I,0 and v I,0 are the position vector and velocity vector respectively;

[0080] Step 4.1.2, click on the arrival point on the map: z impact =[lon impact ,lat impact ,ele impact ] T , where z impact is the geodetic coordinate of the arrival point, lon impact is longitude, lat impact is latitude, ele impact is the elevation;

[0081] Step 4.2, calculate the reachable area of ​​key nodes: record the key nodes including the starting point and the passing points as s I,i , with s I,i is the initial value, and s is calculated by step 3 I,i The reachable area Φ corresponding to the point i , where the arrival point z impact The reachable area Φ at the starting point i Inside;

[0082] Step 4.3, select the next waypoint: in the current reachable area Φ i Select the next waypoint in the i+1 =[lon i+1 ,lat i+1 ,ele i+1 ] T ;

[0083] Step 4.4, calculate the initial value of the route iteration to the next key node, including:

[0084] Step 4.4.1, record the current key node state vector as s I,i The combination of the jth group of turning curvature and turning angle threshold is recorded as K j =[c j ,δ j ] T , with s I,i , Kj As the initial value for extrapolation, calculate the flight route l through step 2 i,j ;

[0085] Step 4.4.2, number the extrapolated initial value parameters:

[0086] K 1,0 =K j , K 3,0 =K j+1 , l 1,0 = l i,j , l 3,0 = l i,j+1 ;

[0087] Step 4.4.3, for K 2,0 Initialization assignment:

[0088] Step 4.4.4, with s I,i , K 2,0 As the initial value, the remaining distance d begins to increase as the extrapolation termination condition, and is brought into step 2 to solve the flight path l 2,0 ;

[0089] Step 4.5, iteratively solving the route to the next key node, including:

[0090] Step 4.5.1, update the extrapolation parameters:

[0091] If i+1 Located in 1,k With l 2,k If the envelope is within the reachable area, the extrapolation parameters are updated as follows:

[0092] K 1,k+1 =K 1,k , K 3,k+1 =K 2,k , l 1,k+1 = l 1,k , l 3,k+1 = l 2,k ;

[0093] If i+1 Located in 2,k With l 3,k If the envelope is within the reachable area, the extrapolation parameters are updated as follows:

[0094] K 1,k+1 =K 2,k , K 3,k+1 =K 3,k , l 1,k+1 = l 2,k , l 3,k+1 = l 2,k ;

[0095] Step 4.5.2, for K 2,k+1 To assign a value:

[0096] Step 4.5.3, with s I,i , K 2,k+1 As the initial value, the remaining distance d starts to increase as the extrapolation termination condition, and is brought into step 2 to solve the flight route l 2,k+1 ;

[0097] Step 4.5.4, repeat steps 4.5.1 to 4.5.3 until z i+1 With route l 2,k+1 The iteration is terminated when the minimum distance of the ground projection is no longer reduced or the minimum distance is less than the set threshold. The current route l 2,k+1 For flight

[0098] The device starts from the key node z i To key node z i+1 flight routes;

[0099] Step 4.5.5, determine z i+1 Is it the arrival point z? impact , such as z i+1 is a waypoint, then let l 2,k+1 The last point is the key node s I,i+1 , bring it into step 4.2 for calculation; if z i+1 To reach point z impact , the calculation terminates and goes to step 4.6;

[0100] Step 4.6, output the entire flight route: connect all the generated flight routes to generate a route from the starting point s of the gliding segment I,0 To the arrival point z impact The complete route l = {l 2,0 ,……,l 2,k+1}.

[0101] Furthermore, step 1.1.1 also includes the following sub-steps:

[0102] Step 1.1.1.1, calculate the state vector of the Earth-centered inertial system of the aircraft: s I =[r I ,v I ] T , where r I is the position vector, r I =[x,y,z] T , velocity vector v I =[v x ,v y ,v z ] T;

[0103] Step 1.1.1.2, calculate the aircraft airspeed vector v a :

[0104]

[0105] Where ω0 is the angular velocity of the Earth's rotation, v wx 、v wy 、v wz are the components of the target local wind speed in the directions of the three coordinate axes of the geocentric inertial system;

[0106] Step 1.1.1.3, calculate the drag acceleration a AD :

[0107]

[0108] In the formula, a HLA It is the acceleration of the aircraft during the gliding flight phase in the atmosphere;

[0109] Step 1.1.1.4, calculate the lift acceleration vector a AL :

[0110] a AL =(a′ AL ·u Ja ) Ja +(a′ AL ·u Ka ) Ka ,

[0111] In the formula, a′ AL is the lift acceleration, a′ AL =a HLA -a G -a AD ,

[0112] Further, step 1.2.2 includes: using the following formula to solve the maximum value c of the mass-to-drag ratio corresponding to the Mach number m: M / Dmax,m :

[0113]

[0114] In the formula, f interp (X,Y,x) is the linear interpolation function of adjacent points, X is the independent variable sequence, Y is the dependent variable sequence, and x is the interpolation point.

[0115] Furthermore, step 2 also includes the following sub-steps:

[0116] Step 2.0.1, obtain the vertical upward direction vector u at the position of the aircraft r :

[0117]

[0118] In the formula, a G,r is the space position of the aircraft in the inertial system at the center of the aircraft I The gravitational acceleration vector at ;

[0119] Step 2.0.2, calculate the direction vector u perpendicular to the ballistic plane j :

[0120]

[0121] Step 2.0.3, calculate the horizontal direction vector u on the ballistic plane i :

[0122]

[0123] Step 2.0.4, calculate the horizontal component v of the inertial system velocity Ii :

[0124] v Ii =(v I ·u i ) i ;

[0125] Step 2.0.5, calculate the angle θ between the inertial system velocity and the horizontal direction I :

[0126]

[0127] Further, step 2.2.1.1 includes: calculating the angle δ between the current velocity vector and the velocity vector at the time of extrapolation initialization by the following formula:

[0128]

[0129] In the formula, v Ii,0 is the horizontal component of the inertial velocity of the aircraft at the initial moment, v Ii,i It is the horizontal component of the inertial system velocity of the aircraft at the current moment.

[0130] Further, step 2.2.1.2 includes: the following formula is used to calculate the centripetal acceleration L required for the aircraft to fly at a curvature C in the horizontal direction: j :

[0131]

[0132] In the formula, v Ii is the horizontal component of the vehicle's inertial velocity.

[0133] Further, step 2.2.3 includes: calculating the lift acceleration L in the equilibrium gliding state by the following formula: k :

[0134]

[0135] Furthermore, in step 2.2.5.3, the lift direction mark D j The calculation is as follows:

[0136]

[0137] Furthermore, step 3.3 includes: taking the aircraft airspeed decreasing to a certain Mach number or the aircraft altitude decreasing to a certain altitude or extrapolating to a certain time as the iteration termination condition.

[0138] Furthermore, step 4.2 also includes: if the arrival point is not located at Φ i If it is within, regenerate the current key node or select a new arrival point.

[0139] The advantages of the present invention are:

[0140] The method for designing the glide segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling of the present invention extracts the maneuvering characteristics of the high-speed gliding aircraft from the non-cooperative measurement data and establishes a characteristic database, performs dynamic modeling on the glide segment of the aircraft, obtains the flight route of the aircraft using the characteristic database and the state vector measured in real time, and then calculates the reachable area of ​​the aircraft, and finally generates the entire glide segment route with the starting point of the glide segment as the initial value and the arrival point as the end point through a number of waypoints. Therefore, the present invention does not need to obtain complex parameters such as the aircraft's angle of attack, roll angle, yaw angle, speed, and rudder angle, as well as the aircraft's real control equations, and can solve the acceleration of the glide process of the high-speed gliding aircraft with only some of the measured aircraft position and speed parameters, and realizes the controllable planning and design of the entire glide segment flight route by autonomously selecting the starting point, arrival point, and waypoints. BRIEF DESCRIPTION OF THE DRAWINGS

[0141] The above and / or additional features and advantages of the present invention will become more easily understood through the following description with reference to the accompanying drawings, in which:

[0142] Figure 1 It is a flow chart of a method for designing a glide segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling of the present invention;

[0143] Figure 2 is the calculation result of the reachable area of ​​the starting point in the example of the present invention;

[0144] Figure 3 is the calculation result of the reachable area of ​​waypoint 1 in the example of the present invention;

[0145] Figure 4 is the calculation result of the reachable area of ​​waypoint 2 in the example of the present invention;

[0146] Figure 5 It is the complete route planning result from the starting point to the arrival point in the example of the present invention;

[0147] Figure 6 It is a comparison of the calculation results of the reachable areas of each key node in the example of the present invention;

[0148] Figure 7 It is the time-elevation curve of the route designed in the example of the present invention;

[0149] Figure 8 It is the time-speed curve of the route designed in the example of the present invention. DETAILED DESCRIPTION

[0150] The present invention will be described in detail below with reference to the accompanying drawings by means of exemplary embodiments of the present invention. It should be noted that the following detailed description of the present invention is only for the purpose of illustration, and is not intended to limit the present invention.

[0151] The present invention provides a method for designing the gliding segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling, which is used to solve the route from the gliding starting point to the waypoint, from the waypoint to the waypoint, and from the waypoint to the arrival point by utilizing a new high-speed gliding aircraft dynamics modeling method, so that only the starting point, the arrival point, and the waypoint need to be selected, and the active route design of the gliding aircraft can be realized without relying on the model design department.

[0152] Reference Figure 1 As an exemplary embodiment of the present invention, a method for designing a glide segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling comprises the following steps:

[0153] Step S1, extracting high-speed gliding aircraft maneuvering characteristics from non-cooperative measurement data and establishing a model maneuvering characteristics database;

[0154] Step S2, establishing a dynamic model for the gliding section of the high-speed gliding aircraft, solving and integrating the acceleration vector using the established model maneuvering characteristic database and the state vector measured in real time, so as to generate the flight route of the aircraft;

[0155] Step S3, designing an aircraft turning curvature sequence and an aircraft turning angle threshold sequence and bringing them into step S2 to obtain multiple flight tracks, thereby obtaining an aircraft reachable area;

[0156] Step S4, taking the starting point of the gliding section as the initial value and the arrival point as the end point, and passing through a number of passing points, generates a full gliding section route.

[0157] Step S1 may include sub-step S1.1: extracting maneuvering features from a single full-flight track data of the high-speed gliding aircraft.

[0158] The key to designing a high-speed glider route is to master its maneuverability. By accumulating and analyzing measured flight data, the maneuverability characteristics can be extracted. The main maneuverability characteristics of an aircraft include drag coefficient, lift-to-drag ratio, and normal overload capacity.

[0159] Acceleration analysis is performed on the smoothed non-cooperative measurement data to extract drag acceleration, lift acceleration and normal overload from a single full-flight track data.

[0160] The acceleration a of a high-speed glider during its gliding flight in the atmosphere HLA It can be decomposed into the gravitational acceleration a G , drag acceleration a AD , lift acceleration a AL sum:

[0161] a HLA =a G +a AD +a AL (1)

[0162] The state vector of the Earth-centered inertial system of the aircraft is expressed as follows, where r I is the position vector, v I is the velocity vector.

[0163]

[0164] Target airspeed vector v a The calculation is as follows:

[0165]

[0166] Where ω0 is the angular velocity of the Earth's rotation, v wx 、v wy 、v wz is the component of the local wind speed in the three coordinate axes of the inertial system. Since the air resistance acceleration is opposite to the airspeed and the lift acceleration is perpendicular to the airspeed acceleration, the drag acceleration can be calculated by the following formula:

[0167]

[0168] After completing the drag acceleration calculation, the lift acceleration can be directly calculated using formula (1) and recorded as a′ AL :

[0169] a′ AL =a HLA -aG -a AD (5)

[0170] However, since the target acceleration data a HLA It is obtained by smoothing the original measurement data and performing second-order difference, so there are some algorithm errors that need to be corrected.

[0171] First calculate the direction vector u perpendicular to the airspeed vector and gravity acceleration vector Ja :

[0172]

[0173] Then calculate the direction vector u which is perpendicular to the airspeed vector and coplanar with the gravity acceleration vector and the airspeed vector Ka :

[0174]

[0175] The corrected lift acceleration vector is expressed as a AL :

[0176] a AL =(a′ AL ·u Ja ) Ja +(a′ AL ·u Ka ) Ka (8)

[0177] g0 is the standard sea level gravity acceleration, and the normal overload is calculated as follows:

[0178]

[0179] The lift-to-drag ratio is calculated as follows:

[0180]

[0181] The mass-to-drag ratio is calculated as follows:

[0182]

[0183] The relationship between the distance d (remaining range) between the aircraft and the target and the aircraft's flight altitude directly reflects the aircraft's control strategy and guidance mode, and is therefore also an important indicator parameter. F is the position vector of the spacecraft's Earth-centered fixed system, r F,Target is the target geocentric fixed position mass, and the remaining range is calculated as follows:

[0184]

[0185] So far, the extraction of all motion features of a single model and a single track has been completed. It can be seen that the extraction of these maneuvering features only depends on some measured aircraft position and velocity parameters, and has no demand for other complex high-order state observations that are difficult to obtain, which is suitable for feature extraction of non-cooperative targets.

[0186] Step S1 also includes sub-step S1.2: establishing a model maneuverability feature database.

[0187] Motion characteristic parameters are accumulated according to the aircraft model, and the mass-to-drag ratio, lift-to-drag ratio, and normal overload characteristics are accumulated respectively. With 1 Mach as the step length, the mean μ and variance σ of the mass-to-drag ratio, lift-to-drag ratio, and normal overload are statistically calculated in each step length. μ+3σ is taken as the maximum value of the characteristic parameter corresponding to the median Mach number of the step length interval, and μ-3σ is taken as the minimum value of the characteristic parameter corresponding to the median Mach number of the step length interval (for example, if μ+3σ=2 is calculated in the interval from 0.5 Mach to 1.5 Mach, it is considered that 1 Mach corresponds to the maximum lift-to-drag ratio of 2). The characteristic parameter maximum and minimum value sequences and the corresponding Mach number sequences are established respectively.

[0188] Establish characteristic parameter sequence, including the maximum value sequence of mass-to-resistance ratio C M / Dmax,Ma , the minimum value sequence of mass-to-resistance ratio C M / D min,Ma , lift-to-drag ratio maximum sequence C L / D max,Ma , lift-to-drag ratio minimum value sequence C L / D min,Ma , Normal overload maximum sequence G Lmax,Ma , Normal overload minimum sequence G Lmin,Ma :

[0189]

[0190] The corresponding Mach number sequence is as follows:

[0191]

[0192] Where m1 is the minimum flight Mach number of this model rounded up, m n The maximum flight Mach number of this model is rounded down.

[0193] Using the adjacent point linear interpolation function, substituting into the characteristic parameter sequence and the Mach number sequence, the maximum value c of the mass-to-drag ratio corresponding to any Mach number m can be solved. M / Dmax,m :

[0194]

[0195] where f interp (X,Y,x) is the linear interpolation function of adjacent points, X is the independent variable sequence, Y is the dependent variable sequence, and x is the interpolation point. Using a method similar to formula (15), the minimum value c of the mass-to-drag ratio corresponding to any Mach number can be completed. M / Dmin,m , lift-to-drag ratio maximum value cL / Dmax,m , lift-to-drag ratio minimum value c L / Dmin,m , Normal overload maximum value g Lmax,m , minimum normal overload g Lmin,m The solution.

[0196] The remaining flight distance and flight altitude data of the same type of aircraft after the downward pressure maneuver at the end of the flight are accumulated. The average altitude value in each step is calculated with the ΔD flight distance as the step length, and it is used as the average flight altitude corresponding to the remaining flight distance at the midpoint of the step length (for example, d i -d i+1 =ΔD, count the remaining distance d i to d i+1 The average flight altitude is h mean,i , then h mean,i As the remaining distance d i +ΔD), and establish the average flight height sequence at the end of the flight:

[0197] H mean,D =[h mean,1 ,…,h mean,n ] (16)

[0198] The corresponding remaining flight sequence is as follows:

[0199]

[0200] By using the linear interpolation function of adjacent points and substituting into the remaining flight sequence and the average flight altitude sequence, we can use formula (18) to solve the average flight altitude h corresponding to any remaining flight distance d: mean,d , where f interp is the linear interpolation function of adjacent points, and its usage is introduced above.

[0201]

[0202] For step S2, the most important motion state of a high-speed gliding aircraft in the gliding stage is balanced gliding, so it is necessary to perform dynamic modeling of the balanced gliding state. Step S2 includes sub-steps S2.1: solving the angle θ between the aircraft speed and the horizontal direction a .

[0203] The vertical upward direction vector u of the aircraft's location r It is expressed as formula (19), where a G,r is the space position of the aircraft in the inertial system at the center of the aircraft I The gravitational acceleration vector at .

[0204]

[0205] Due to the influence of non-spherical gravity, the center of the aircraft's flight is not the center of the earth. Under the second-order gravity model, the center vector R g The calculation method is shown in formula (20):

[0206]

[0207] Direction vector u perpendicular to the flight plane j The calculation is as follows:

[0208]

[0209] The horizontal direction vector u on the flight path plane i The calculation is as follows:

[0210]

[0211] The direction vector u is on the ballistic plane and perpendicular to the velocity direction k The calculation is as follows:

[0212]

[0213] Horizontal component of inertial system velocity v Ii The calculation is as follows:

[0214] v Ii =(v I ·u i ) i (twenty four)

[0215] The angle θ between the inertial system velocity and the horizontal direction I The calculation is shown in formula (25). When the velocity pitch angle is upward, θ I The value is positive, the speed pitch angle is downward θ I The value is negative.

[0216]

[0217] The angle between airspeed and horizontal direction θ a The calculation is shown in formula (26):

[0218]

[0219] Step S2 also includes sub-step S2.2: establishing a dynamic model and solving the acceleration vector.

[0220] Specifically, step S2.2 includes sub-step S2.2.1: determining the motion state of the aircraft.

[0221] The component a of the acceleration due to gravity on the plane where the lift vector of the aircraft lies GKJ :

[0222] a GKJ =||a G ||cosθ a (27)

[0223] Use equation (28) to calculate the angle δ between the current velocity vector and the extrapolated initial velocity vector, where v Ii,0 is the horizontal component of the inertial velocity of the aircraft at the initial moment, v Ii,i is the horizontal component of the inertial velocity of the aircraft at the current moment, and its calculation method is shown in formula (24).

[0224]

[0225] Use formula (29) to calculate the centripetal acceleration L required for the aircraft to fly in the horizontal direction with a curvature C: j Where δ threshold is the turning angle threshold. When the turning angle of the aircraft relative to the extrapolated initial state reaches δ threshold Then it starts to fly level.

[0226]

[0227] Design lift direction mark D j :

[0228]

[0229] Set the angle threshold between the aircraft inertial system velocity and the horizontal direction in the balanced gliding state to θ ST , and design three discriminants as follows:

[0230] θ I <θ ST (31)

[0231] θ ST <θ I ≤0 (32)

[0232] θ I >0 (33)

[0233] When condition (31) is met, it is determined that the aircraft is in the downward pressure phase of the trajectory and is about to perform a pull-up operation, and the process proceeds to step S2.2.2 to be described below; when condition (32) is met, it is determined that the aircraft is in a balanced gliding state and the process proceeds to step S2.2.3 to be described below; when condition (33) is met, it is determined that the aircraft is in a ballistic flight state before entering a balanced gliding state and the process proceeds to step S2.2.4 to be described below.

[0234] Step S2.2 includes sub-step S2.2.2: calculating the lift acceleration in the pull-up state.

[0235] First, determine whether its current elevation is lower than the elevation of the pull-up starting point (the elevation of the pull-up starting point is obtained by counting the actual measured data of the corresponding model).

[0236] If the aircraft altitude is higher than the pull-up point altitude, then the instantaneous ballistic surface lift acceleration L k =0; if the aircraft altitude is less than or equal to the pull-up point altitude, then L k Calculated by formula (34); where a GKJ Calculated by formula (27); g Lmax,m is the maximum overload of the aircraft at the current Mach number, and its calculation method is shown in step S1.

[0237] L k =a GKJ +g Lmax,m ·g0 (34)

[0238] Step S2.2 also includes sub-step S2.2.3: calculating the lift acceleration in the equilibrium gliding state.

[0239] The magnitude of the instantaneous ballistic surface lift acceleration L k Calculated by formula (35):

[0240]

[0241] In addition, step S2.2 also includes sub-step S2.2.4: calculating the ballistic flight lift acceleration.

[0242] At this time, there is no lift acceleration, and the L k , L j To assign a value:

[0243]

[0244] Step S2.2 also includes sub-step S2.2.5: determining whether flight altitude control is required.

[0245] Use formula (12) to calculate the remaining distance d from the aircraft to the target and the current flight altitude h t , substitute the remaining distance d into equation (18) to solve the corresponding average flight altitude h mean,d If the flight altitude control is not started and the current h t <h mean,d Then the flight altitude control is not performed and the process goes to step S2.6; otherwise, the flight altitude control is started.

[0246] Use formula (38) to calculate the lift acceleration L of the ballistic surface k , and update the coefficient h int,i The value of c A 、cB 、c c is the adjustment coefficient with a fixed value in advance; θ I is the angle between the inertial system velocity and the horizontal direction, calculated by formula (25); h i-1 h is the altitude of the aircraft during the last iteration calculation; int,i Initialize h for the height control integral correction int,0 =0, and then update its value using equation (38) in each iteration.

[0247]

[0248] Step S2.2 also includes sub-step S2.2.6: synthesizing and correcting the lift acceleration to a total lift acceleration L. After respectively calculating the lift acceleration of the pull-up state in step S2.2.2, the lift acceleration of the balanced gliding state in step S2.2.3, and the lift acceleration of the ballistic flight in step S2.2.4, step S2.2.5 is executed. If it is determined that no altitude control is required, the process directly jumps to step S2.2.6, otherwise the flight altitude control in step S2.2.5 is executed.

[0249] In step S2.2.6, the maximum lift acceleration a under the current Mach number can be obtained by using equation (39) in the current state of the aircraft: ALmax,m , where c M / Dmin,m is the minimum mass-to-drag ratio at the current Mach number, c L / Dmax,m is the maximum lift-to-drag ratio at the current Mach number, c M / Dmin,m 、c L / Dmax,m It is obtained by querying the model maneuverability characteristics database, and its calculation method is shown in step S1.2.

[0250]

[0251] If L k >a ALmax,m , the total lift that can be provided at this time is less than the expected ballistic surface lift, then:

[0252]

[0253] If L k ≤a ALmax,m , at this time, the pull-up lift that can be provided is greater than the expected pull-up lift, and the remaining lift can be used for horizontal maneuvers. At this time, the maximum lift acceleration that can be used for horizontal maneuvers is:

[0254]

[0255] If the current |L j |>L jmax , then let L j =Dj L jmax , otherwise L j No change.

[0256] Calculate the resultant force L on the projectile perpendicular to the airspeed direction ki :

[0257]

[0258] Calculate the maximum normal overload g at the current Mach number through step S1 Lmax,m , if L ki >g Lmax,m , then recalculate L j :

[0259]

[0260] Calculate the total lift acceleration L:

[0261]

[0262] Aircraft roll angle The calculation is as follows:

[0263]

[0264] Step S2.2 also includes sub-step S2.2.7: Calculate the drag acceleration a AD .

[0265] In the shorter Mach number range, it is assumed that the lift and drag acceleration have linear characteristics. Use equation (39) to calculate the maximum drag acceleration a at the current Mach number: ADmax,m and lift acceleration a ADmax,m , and then calculate the minimum drag acceleration a at the current Mach number through formula (46) ADmin,m and lift acceleration a ADmin,m 、c M / Dmax,m 、c L / Dmin,m The calculation method of is shown in step S1.2.

[0266]

[0267] Design the following equation to complete the resistance acceleration a AD Calculation:

[0268]

[0269] Step S2.2 also includes sub-step S2.2.8: acceleration vector synthesis and integration.

[0270] After completing the calculation of the lift acceleration and drag acceleration, calculate the drag and lift acceleration vectors of the aircraft at this moment:

[0271]

[0272] Combining formula (1) and (2), we can know that:

[0273]

[0274] right The state vector s can be realized by integration I If the iteration termination condition is met, the iteration is terminated, otherwise, the calculation is continued by jumping to step S2.2.1.

[0275] Step S3 can realize the prediction of the reachable range of the aircraft. The specific calculation process of step S3 is as follows:

[0276] Design curvature sequence C:

[0277] C=[c1,…,c n＝2m+3 ],m∈N + (50)

[0278]

[0279] Where, d cr,j is the intermediate variable function, the subscript j is the independent variable of the function, c r2 is the maximum turning radius, c r1 is the minimum turning radius.

[0280] Design turning angle threshold sequence δ threshold :

[0281]

[0282] The turning angle threshold and the turning curvature are constructed into a state vector combination K:

[0283] K i =[c i ,δ i ] T ,i=1…n (52)

[0284] K i Substitute into step S2 in sequence, taking the state vector s at a certain moment I,t As the initial point of extrapolation, the iteration termination condition is when the aircraft airspeed drops to a certain Mach number or the aircraft altitude drops to a certain altitude or when the aircraft is extrapolated to a certain time. The n flight routes generated under each combination of curvature and turning angle thresholds are calculated, and the adjacent route endpoints and the two outermost flight routes are connected in sequence to form a closed envelope area. The envelope ground projection area is the reachable area of ​​the aircraft, which is recorded as Φ t .

[0285] In step S4, the glide segment starting point is used as the initial value, the arrival point is used as the end point, and a number of passing points are passed to generate the entire glide segment route. The starting point, passing points, and arrival points can be collectively referred to as key nodes, and the key nodes need to be constructed first.

[0286] Step S4 includes step S4.1: constructing a starting point and an arrival point.

[0287] As shown in formula (53), the state vector of the inertial system at the initial point of the aircraft gliding phase is denoted as s I,0 , where the second digit of the superscript is the point number.

[0288] s I,0 =[r I,0 ,v I,0 ] T (53)

[0289] s I,0 The measured data can be used directly, or generated by translating and rotating the measured data.

[0290] The arrival point is selected on the map according to actual needs, as shown in formula (54), and the geodetic coordinates of the arrival point are marked as z impact , from the longitude lon impact , latitude lat impact 、Elevation impact constitute.

[0291] z impact =[lon impact ,lat impact ,ele impact ] T (54)

[0292] Step S4 also includes step S4.2: calculating the reachable area of ​​the key node.

[0293] The current key node, whether it is the starting point or the passing point, is recorded as s I,i , with s I,i is the initial value, and s is calculated by the method in step S3. I,i The reachable area Φ corresponding to the point i Arrival point z impact Must be located in the reachable area Φ of the starting point i Inside, and satisfy the following formula:

[0294] [lon impact ,lat impact ] T ∈Φ i (55)

[0295] If the arrival point is not located at Φ iIf the current key node s is within , the next step of calculation cannot be performed and the current key node s needs to be regenerated. I,i , or select a new arrival point z impact , and is calculated through step S4.2 until the condition (55) is met.

[0296] To directly solve the route to the arrival point, go to step S4.4; to solve the route to the next waypoint, go to step S4.3 which will be described below.

[0297] Step S4 also includes step S4.3: selecting the next waypoint.

[0298] In the currently accessible area Φ i Select the next waypoint in the list. The geodetic coordinates of this point are marked as z i+1 =[lon i+1 ,lat i+1 ,ele i+1 ] T , and this point satisfies the following formula:

[0299] [lon i+1 ,lat i+1 ] T ∈Φ i (56)

[0300] Step S4 also includes step S4.4: calculating the initial value of the route iteration to the next key node.

[0301] The current key node state vector is s I,i , the combination of the jth group of turning curvature and turning angle threshold is recorded as K j =[c j ,δ j ] T ; I,i , K j The flight path calculated in step S2 is used as the initial value for extrapolation and is denoted as l i,j During the calculation process, if the next key node is a waypoint, no altitude control is performed; if the next key node is an arrival point, altitude control is performed. i,j , l i,j+1 The enclosed reachable area is denoted as Φ i,j .

[0302] The geodetic coordinates of the next key node are marked as z i+1 =[lon i+1 ,lat i+1 ,ele i+1 ] T , if formula (57) is satisfied, the flight route of the aircraft from the current key node to the next key node can be calculated.

[0303] [loni+1 ,lat i+1 ] T ∈Φ i,j (57)

[0304] To facilitate calculation, the extrapolated initial value parameters are renumbered using formula (58), K 1,0 , K 3,0 , l 1,0 , l 3,0 The first digit of several parameter subscripts is the parameter number, and the second digit is the number of iterations.

[0305]

[0306] K 2,0 Initialization assignment is as follows:

[0307]

[0308] In s I,i , K 2,0 As the initial value, the remaining distance d starts to increase as the extrapolation termination condition (d is calculated by formula (12)), and is brought into step S2 to solve the flight route l 2,0 During the calculation process, if the next key node is a waypoint, no altitude control is performed; however, if the next key node is an arrival point, altitude control is performed.

[0309] Step S4 also includes step S4.5: iteratively solving the route to the next key node.

[0310] In the kth iteration, z i+1 If it is located at l 1,k With l 2,k If the envelope is within the reachable area, the extrapolation parameters are updated as follows:

[0311]

[0312] z i+1 If it is located at l 2,k With l 3,k If the envelope is within the reachable area, the extrapolation parameters are updated as follows:

[0313]

[0314] K 2,k+1 To assign a value:

[0315]

[0316] In s I,i , K 2,k+1 As the initial value, the remaining distance d starts to increase as the extrapolation termination condition, and is brought into step S2 to solve the flight route l2,k+1 .

[0317] Repeat step S4.5 until z i+1 With route l 2,k+1 If the minimum distance of the ground projection of is no longer reduced or is less than a certain threshold, the iteration is terminated. 2,k+1 The route is the aircraft from the key node z i To key node z i+1 flight routes.

[0318] Such as i+1 is a waypoint, then let l 2,k+1 The last point is the key node s I,i+1 , bring it into step S4.2 for calculation.

[0319] Such as i+1 To reach point z impact , the calculation terminates and goes to the following step S4.6.

[0320] Step S4 also includes step S4.6: outputting the entire flight route.

[0321] Connecting the generated routes of all segments will generate a route from the starting point s I,0 To the arrival point z impact The complete route:

[0322] l={l 2,0 ,……,l 2,k+1} (63)

[0323] As described above, the method for designing the glide segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling of the present invention extracts the maneuvering characteristics of the high-speed gliding aircraft from the non-cooperative measurement data and establishes a characteristic database, and also performs dynamic modeling on the glide segment of the aircraft, obtains the flight route of the aircraft using the characteristic database and the state vector measured in real time, and then calculates the reachable area of ​​the aircraft, and finally generates the entire glide segment route with the starting point of the glide segment as the initial value and the arrival point as the end point through a number of waypoints. Therefore, the present invention does not need to obtain complex parameters such as the aircraft's angle of attack, roll angle, yaw angle, speed, and rudder angle, as well as the aircraft's real control equations, and can solve the acceleration of the glide process of the high-speed gliding aircraft with only some of the measured aircraft position and speed parameters, and realizes the controllable planning and design of the entire glide segment flight route by autonomously selecting the starting point, arrival point, and waypoints.

[0324] Next, refer to Figures 2 to 8 The method for designing the glide segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling provided by the present invention is further described with reference to an example.

[0325] In this example, firstly, a starting state vector of the gliding segment of a high-speed gliding aircraft is constructed and generated based on reasonable assumptions as the starting point of the route design, and a point on the map is selected as the arrival point.

[0326] Then, use step S4.2 to calculate the reachable area of ​​the aircraft at the starting point. When calculating the reachable area, bring it into the calculation process of step S3, and c in formula (50) r1 The value is 100km, c r2 The value is 10000km, m is 7, n is km17, and the iteration termination condition is when the aircraft Mach number drops below Mach 2. Figure 2 The calculation result of the reachable area of ​​the aircraft at the starting point is shown. The dotted area in the figure is the reachable area, ◇ is the starting point, and ▽ is the arrival point. It can be seen that the arrival point is within the reachable area of ​​the starting point, and the route design can continue.

[0327] Then, use step S4.3 to select waypoint 1 in the current reachable area, complete the route calculation from the starting point to waypoint 1 through steps S4.4 and S4.5, and then use step S4.2 to calculate the reachable area when the aircraft reaches waypoint 1. Figure 3 This is the calculation result of the reachable area of ​​the aircraft at waypoint 1. The ○ in the figure indicates the position of the waypoint. It can be seen that the arrival point is within the reachable area, and the route design can continue.

[0328] Then, step S4.3 is used to select waypoint 2 in the current reachable area, and the route calculation from waypoint 1 to waypoint 2 is completed through steps S4.4 and S4.5, and then step S4.2 is used to calculate the reachable area when the aircraft reaches waypoint 2. Figure 4 The calculation result of the reachable area of ​​the aircraft at waypoint 2 is shown. It can be seen from the figure that the arrival point is still within the reachable area and the route design can continue.

[0329] Then, the arrival point is taken as the target key point, and the route calculation from the waypoint 2 to the arrival point is completed through steps S4.4 and S4.5. The parameter values ​​used in formula (38) are as follows:

[0330]

[0331] like Figure 5 The figure shows the complete designed route of the aircraft from the starting point to the arrival point. It can be seen that the route has a complex shape and passes through waypoints 1 and 2, showing that the aircraft has strong maneuverability.

[0332] Comparison of the calculation results of the reachable area when the aircraft reaches each key node Figure 6 It can be seen that the calculation result of the reachable area of ​​each key node is a true subset of the calculation result of the reachable area of ​​the previous key node.

[0333] Figure 7 and Figure 8 They are the time-altitude curve and time-speed curve of the complete designed route. It can be observed from these two figures that the process from the starting point to the waypoint 2 is mainly a balanced gliding state, while the process from the waypoint 2 to the arrival point is actively controlled in terms of altitude using the relevant parameters of the model maneuvering characteristics database.

[0334] In summary, this example verifies the effectiveness of the high-speed gliding aircraft gliding segment route design method based on non-cooperative measurement data modeling provided by the present invention.

[0335] Finally, it should be noted that the features mentioned in the above description of the exemplary embodiments can be combined with other embodiments, combined with the features in other embodiments or replace the corresponding features in other implementations. The technical solutions obtained by these combinations or replacements should also be deemed to be included in the protection scope of the present invention.

Claims

1. A method for designing a glide segment route for a high-speed gliding aircraft based on non-cooperative measurement data modeling, characterized in that: The following steps are involved: Step 1, extracting high-speed gliding aircraft maneuvering characteristics from non-cooperative measurement data and establishing a model maneuvering characteristics database, includes the following sub-steps: Step 1.1, extract maneuvering features from a single full-flight track data of a high-speed glider aircraft, the maneuvering features include the mass-to-drag ratio g M / D , lift-to-drag ratio g L / D , normal overload g L , include: Step 1.1.1, calculate the normal overload g L : Where g0 is the standard sea level gravity acceleration, a AL is the lift acceleration vector; Step 1.1.2, calculate lift-to-drag ratio g L / D : In the formula, a AD is the drag acceleration; Step 1.1.3, calculate the mass-to-resistance ratio g M / D : In the formula, ρ r is the atmospheric density at the location of the aircraft, v a is the aircraft airspeed vector; Step 1.1.4, calculate the distance d between the aircraft and the target: In the formula, r F is the position vector of the Earth-centered fixed system of the aircraft, r F,Target is the target geocentric fixed position mass; Step 1.2, establish a model maneuverability feature database, including: Step 1.2.1, establish a characteristic parameter sequence, including: Mass-to-resistance ratio maximum sequence C M / Dmax,Ma =[c M / Dmax,m1 ,…,c M / Dmax,mn ]、 Mass-to-resistance ratio minimum sequence C M / Dmin,Ma =[c M / Dmin,m1 ,…,c M / Dmin,mn ]、 Lift-to-drag ratio maximum sequence C L / Dmax,Ma =[c L / Dmax,m1 ,…,c L / Dmax,mn ]、 Lift-to-drag ratio minimum value sequence C L / Dmin,Ma =[c L / Dmin,m1 ,…,c L / Dmin,mn ]、 Normal overload maximum sequence G Lmax,Ma =[g Lmax,m1 ,…,g Lmax,mn ]、 Normal overload minimum sequence G Lmin,Ma =[g Lmin,m1 ,…,g Lmin,mn ]、 And the average flight altitude sequence H at the end of the flight mean,D =[h mean,d1 ,…,h mean,dn ], Where Ma={m1,…,m n } is the Mach number sequence, D={d1,…,d n } is the remaining flight sequence; Step 1.2.2, use the linear interpolation method to solve the maximum value c of the mass-to-drag ratio corresponding to the Mach number m M / Dmax,m , minimum value of mass-to-resistance ratio c M / Dmin,m , lift-to-drag ratio maximum value c L / Dmax,m , lift-to-drag ratio minimum value c L / Dmin,m , Normal overload maximum value g Lmax,m , minimum normal overload g Lmin,m , and the average flight altitude h corresponding to the remaining flight distance d mean,d ; Step 2, establishing a dynamic model for the glide segment of the high-speed glider, using the established model maneuvering characteristic database and the real-time measured state vector to solve and integrate the acceleration vector to generate the flight route of the aircraft, including the following sub-steps: Step 2.1: Solve for the angle θ between the aircraft speed and the horizontal direction a : In the formula, u r is the vertical upward direction vector of the aircraft's location; Step 2.2, establishing a dynamic model and solving the acceleration vector, includes the following steps: Step 2.2.1, aircraft motion state determination, including: Step 2.2.1.1, calculate the angle δ between the current velocity vector and the velocity vector at the time of extrapolation initialization; Step 2.2.1.2, calculate the centripetal acceleration L required for the aircraft to fly at a curvature C in the horizontal direction j ; Step 2.2.1.3, determine the motion state of the aircraft, when θ I <θ ST When θ ST <θ I ≤0, the aircraft is judged to be in a balanced gliding state. I >0, the aircraft is judged to be in a ballistic flight state before entering a balanced gliding state, where θ ST is the angle threshold between the inertial system velocity and the horizontal direction of the aircraft in the balanced gliding state, θ I is the angle between the inertial system velocity and the horizontal direction of the aircraft at the current moment; Step 2.2.2, calculate the lift acceleration L in the pull-up state k : L k =a GKJ +g Lmax,m ·g0, In the formula, a GKJ is the component of gravity acceleration on the plane where the lift vector of the aircraft is located, a GKJ =||a G ||cosθ a , a G is the acceleration due to gravity, g Lmax,m is the maximum overload of the aircraft at the current Mach number; Step 2.2.3: Calculate the lift acceleration L in the equilibrium gliding state k ; Step 2.2.4, calculate the lift acceleration L of ballistic flight k ; Step 2.2.5, determine whether flight altitude control is required: calculate the remaining distance d between the aircraft and the target, and calculate the current flight altitude h t , solve the corresponding average flight height h mean,d , if the current flight altitude is h t <h mean,d , then no flight altitude control is performed and go to step 2.2.6; otherwise, the following flight altitude control is performed: Calculate the lift acceleration L on the trajectory surface k , and update the elevation control integral correction h int,i Values: In the formula, c A 、c B 、c c h is the adjustment coefficient with a fixed value in advance. i-1 is the aircraft altitude at the last iteration, h i is the current altitude of the aircraft; Step 2.2.6, synthesizing and correcting the lift acceleration to the total lift acceleration L, includes: Step 2.2.6.1, calculate the maximum lift acceleration a that can be obtained at the current Mach number m of the aircraft in its current state ALmax,m : a ALmax,m =c L / Dmax,m ·a ADmax,m , In the formula, c M / Dmin,m is the minimum mass-to-drag ratio at the current Mach number, c L / Dmax,m is the maximum lift-to-drag ratio at the current Mach number; Step 2.2.6.2, if L k >a ALmax,m , then let: L k =a ALmax,m , L j =0; if L k ≤a ALmax,m , then calculate the maximum lift acceleration L that can be used for horizontal maneuvers jmax : Step 2.2.6.3, determine the current lift acceleration Lj, if |L j |>L jmax , then let L j =D j L jmax , D j is the lift direction sign, otherwise L j No change; Step 2.2.6.4, calculate the resultant force L on the projectile perpendicular to the airspeed direction ki : Step 2.2.6.5, in L ki >g Lmax,m When L k , L j : Step 2.2.6.6, calculate the total lift acceleration L: Step 2.2.7, calculate the drag acceleration a AD : Step 2.2.7.1, calculate the minimum drag acceleration a at the current Mach number ADmin,m and lift acceleration a ADmin,m : a ALmin,m =c L / Dmin,m ·a ADmin,m , In the formula, c M / Dmax,m is the maximum mass-to-drag ratio at the current Mach number, c L / Dmin,m is the minimum lift-to-drag ratio at the current Mach number; Step 2.2.7.2, calculate the drag acceleration a AD : Step 2.2.8, acceleration vector synthesis and integration, includes: Step 2.2.8.1, calculate the aircraft drag acceleration vector a AD and the lift acceleration vector a AL : a AL =L k ·and ka +L j ·and ja , In the formula, u Ka is the direction vector perpendicular to the airspeed vector and coplanar with the gravity acceleration vector and the airspeed vector, u Ja is the direction vector perpendicular to the airspeed vector and gravity acceleration vector; Step 2.2.8.2, for Integrate and generate the aircraft trajectory, and terminate the iteration when the iteration termination condition is met: In the formula, v I is the velocity vector; Step 3, designing the aircraft turning curvature sequence and the aircraft turning angle threshold sequence and bringing them into step 2 to obtain multiple flight tracks, thereby obtaining the aircraft reachable area, including the following sub-steps: Step 3.1, design the curvature sequence C of the aircraft: C=[c1,…,c n＝2m+3 ],m∈N + , Where, d cr,j is the intermediate variable function, the subscript j is the independent variable of the function, c r2 is the maximum turning radius, c r1 is the minimum turning radius; Step 3.2, design the aircraft turning angle threshold sequence δ threshold : d threshold =[δ1,…,δ n＝2m+3 ],m∈N + , Step 3.3: Combine the turning curvature and the turning angle threshold K i =[c i ,δ i ] T , i=1…n are successively brought into step 2, and the state vector s at a certain moment is I,t As the initial point of extrapolation, perform integral calculation on each flight route, terminate the iteration when the iteration termination condition is met, and continue until all routes are calculated, thereby generating n flight routes, and sequentially connecting the endpoints of adjacent routes with the two outermost flight routes to form a closed envelope area, the ground projection area of ​​which is the reachable area of ​​the aircraft; Step 4, taking the starting point of the gliding segment as the initial value and the arrival point as the end point, and passing through a number of passing points, generates the entire gliding segment route, including the following sub-steps: Step 4.1, construct the starting point and the arrival point, including: Step 4.1.1, select the starting point according to the measured data: s I,0 =[r I,0 ,v I,0 ] T , where s I,0 is the inertial system state vector at the start point of the aircraft glide segment, r I,0 and v I,0 are the position vector and velocity vector respectively; Step 4.1.2, click on the arrival point on the map: z impact =[lon impact ,lat impact ,ele impact ] T , where z impact is the geodetic coordinate of the arrival point, lon impact is longitude, lat impact is latitude, ele impact is the elevation; Step 4.2, calculate the reachable area of ​​key nodes: record the key nodes including the starting point and the passing points as s I,i , with s I,i is the initial value, and s is calculated by step 3 I,i The reachable area Φ corresponding to the point i , where the arrival point z impact The reachable area Φ at the starting point i Inside; Step 4.3, select the next waypoint: in the current reachable area Φ i Select the next waypoint in the i+1 =[lon i+1 ,lat i+1 ,ele i+1 ] T ; Step 4.4, calculate the initial value of the route iteration to the next key node, including: Step 4.4.1, record the current key node state vector as s I,i The combination of the jth group of turning curvature and turning angle threshold is recorded as K j =[c j ,δ j ] T , with s I,i , K j As the initial value for extrapolation, calculate the flight route l through step 2 i,j ; Step 4.4.2, number the extrapolated initial value parameters: K 1,0 =K j ,K 3,0 =K j+1 ,l 1,0 =l i,j ,l 3,0 =l i,j+1 ; Step 4.4.3, for K 2,0 Initialization assignment: Step 4.4.4, with s I,i , K 2,0 As the initial value, the remaining distance d begins to increase as the extrapolation termination condition, and is brought into step 2 to solve the flight path l 2,0 ; Step 4.5, iteratively solving the route to the next key node, including: Step 4.5.1, update the extrapolation parameters: If i+1 Located in 1,k With l 2,k If the envelope is within the reachable area, the extrapolation parameters are updated as follows: K 1,k+1 =K 1,k ,K 3,k+1 =K 2,k ,l 1,k+1 =l 1,k ,l 3,k+1 =l 2,k ; If i+1 Located in 2,k With l 3,k If the envelope is within the reachable area, the extrapolation parameters are updated as follows: K 1,k+1 =K 2,k ,K 3,k+1 =K 3,k ,l 1,k+1 =l 2,k ,l 3,k+1 =l 2,k ; Step 4.5.2, for K 2,k+1 To assign a value: Step 4.5.3, with s I,i , K 2,k+1 As the initial value, the remaining distance d starts to increase as the extrapolation termination condition, and is brought into step 2 to solve the flight route l 2,k+1 ; Step 4.5.4, repeat steps 4.5.1 to 4.5.3 until z i+1 With route l 2,k+1 The iteration is terminated when the minimum distance of the ground projection is no longer reduced or the minimum distance is less than the set threshold. The current route l 2,k+1 For the aircraft from the key node z i To key node z i+1 flight routes; Step 4.5.5, determine z i+1 Is it the arrival point z? impact , such as z i+1 is a waypoint, then let l 2,k+1 The last point is the key node s I,i+1 , bring it into step 4.2 for calculation; if z i+1 To reach point z impact , the calculation is terminated and goes to step 4.6; Step 4.6, output the entire flight route: connect all the generated flight routes to generate a route from the starting point s of the gliding segment I,0 To the arrival point z impact The complete route l = {l 2,0 ,……,l 2,k+1 }.

2. The method for designing a glide segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling according to claim 1, characterized in that: The step 1.1.1 also includes the following sub-steps: Step 1.1.1.1, calculate the state vector of the Earth-centered inertial system of the aircraft: s I =[r I ,v I ] T , where r I is the position vector, r I =[x,y,z] T , velocity vector v I =[v x ,v y ,v z ] T ; Step 1.1.1.2, calculate the aircraft airspeed vector v a : Where ω0 is the angular velocity of the Earth's rotation, v wx 、v wy 、v wz are the components of the target local wind speed in the directions of the three coordinate axes of the geocentric inertial system; Step 1.1.1.3, calculate the drag acceleration a AD : In the formula, a HLA It is the acceleration of the aircraft during the gliding flight phase in the atmosphere; Step 1.1.1.4, calculate the lift acceleration vector a AL : a AL =(a′ AL ·u Ja )u Ja +(a′ AL ·u Ka )u Ka , In the formula, a′ AL is the lift acceleration, a′ AL =a HLA -a G -a AD , 3. The method for designing a glide segment route for a high-speed gliding aircraft based on non-cooperative measurement data modeling according to claim 1 or 2, characterized in that: The step 1.2.2 includes: using the following formula to solve the maximum value c of the mass-to-drag ratio corresponding to the Mach number m M / Dmax,m : In the formula, f interp (X,Y,x) is the linear interpolation function of adjacent points, X is the independent variable sequence, Y is the dependent variable sequence, and x is the interpolation point.

4. The method for designing a glide segment route for a high-speed gliding aircraft based on non-cooperative measurement data modeling according to claim 1 or 2, characterized in that: The step 2 also includes the following sub-steps: Step 2.0.1, obtain the vertical upward direction vector u at the position of the aircraft r : In the formula, a G,r is the space position of the aircraft in the inertial system at the center of the aircraft I The gravitational acceleration vector at ; Step 2.0.2, calculate the direction vector u perpendicular to the ballistic plane j : Step 2.0.3, calculate the horizontal direction vector u on the ballistic plane i : Step 2.0.4, calculate the horizontal component v of the inertial system velocity Ii : v Ii =(v I ·u i )u i ; Step 2.0.5, calculate the angle θ between the inertial system velocity and the horizontal direction I :

5. The method for designing a glide segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling according to claim 1 or 2, characterized in that: The step 2.2.1.1 includes: calculating the angle δ between the current velocity vector and the velocity vector at the initial time of extrapolation by the following formula: In the formula, v Ii,0 is the horizontal component of the inertial velocity of the aircraft at the initial moment, v Ii,i It is the horizontal component of the inertial system velocity of the aircraft at the current moment.

6. The method for designing a glide segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling according to claim 1 or 2, characterized in that: The step 2.2.1.2 includes: the following formula is used to calculate the centripetal acceleration L required for the aircraft to fly at a curvature C in the horizontal direction: j : In the formula, v Ii is the horizontal component of the vehicle's inertial velocity.

7. The method for designing a glide segment route for a high-speed gliding aircraft based on non-cooperative measurement data modeling according to claim 1 or 2, characterized in that: The step 2.2.3 includes: calculating the lift acceleration L in the equilibrium gliding state by the following formula k :

8. The method for designing a glide segment route for a high-speed gliding aircraft based on non-cooperative measurement data modeling according to claim 1 or 2, characterized in that: In step 2.2.5.3, the lift direction mark D j The calculation is as follows:

9. The method for designing a glide segment route for a high-speed gliding aircraft based on non-cooperative measurement data modeling according to claim 1 or 2, characterized in that: The step 3.3 includes: taking the aircraft airspeed decreasing to a certain Mach number or the aircraft altitude decreasing to a certain altitude or extrapolating to a certain time as the iteration termination condition.

10. The method for designing a glide segment route of a high-speed gliding aircraft based on non-cooperative measurement data modeling according to claim 1 or 2, characterized in that: The step 4.2 further includes: if the arrival point is not located at Φ i If it is within, regenerate the current key node or select a new arrival point.

Images