Simulation Platform for SINS/GPS Integrated Navigation System of Hypersonic Vehicles Based on Flight Mechanics

In this study, a simulation platform for an integrated navigation algorithm for hypersonic vehicles based on flight mechanics is designed. In addition, the generation method of inertial measurement unit data and satellite receiver data is introduced. First, the interface relationship between a high-precision six-degree-of-freedom (6DoF) model and the simulation platform in the launch-centered Earth-fixed frame is introduced. Three-axis theoretical specific force and angular velocity are output by the 6DoF model. Accelerometer and gyroscope error models are added, and integral processing of the specific force and angular velocity is performed to obtain velocity increment of the accelerometer and the angular increment of the gyroscope. These data are quantified to obtain the accelerometer and gyroscope pulses. The satellite’s pseudo-range and pseudo-range rate as well as its position and velocity are obtained from the theoretical position, velocity, the attitude of the hypersonic vehicle’s 6DoF model output, and the global positioning system (GPS) satellite broadcast ephemeris. The simulation data can be used for the verification of the loose and tight coupling integrated navigation algorithms. The simulation test verifies the accuracy of the designed method.


Introduction
Near space is recognized as the altitude ranging 20-100 km, which is high for airplanes and too low for satellites. A hypersonic vehicle can travel Mach 5 (five times the speed of sound) or higher in near space [1,2]. Hypersonic vehicles employ strapdown inertial navigation systems (SINS) to obtain comprehensive navigation information, high autonomy, and a high update rate. The X-43A hypersonic vehicle employs an integrated navigation system with a SINS and a global positioning system (GPS) [3]. The X-51A navigation system is equipped with an inertial measurement unit (IMU) and a GPS receiver [4]. The Hypersonic Technology Vehicle 2 operates with a tightly coupled GPS/IMU system for guidance and approximately 3 m of error in the endo-atmospheric glide [5]. The navigation system of the German SHEFEX-2 hypersonic vehicle integrates an IMU, a GPS receiver, and a star tracker [6]. G. Hu et al present a new innovation orthogonality-based robust unscented Kalman filter (IO-RUKF) for INS/GNSS tightly coupled integrated navigation of hypersonic aircraft [7]. The Hypersonic Conventional Strike Weapon utilizes a GPS/INS system for navigation and terminal guidance. Countries such as Russia, Japan, and India have also installed SINS/global navigation satellite system (GNSS) integrated navigation systems in their hypersonic vehicles [1].

Integrated Navigation Simulation Platform
Semi-physical simulation of aircraft is conducted using the HWIL simulation method, which is one of the important aspects in the development of aircraft. The flight control system of a hypersonic vehicle is connected with the HWIL simulation system to reproduce the flight environment of the hypersonic vehicle in air as realistically as possible in underground laboratory conditions and to verify and evaluate the performance indicators of the flight control system.
The overall architecture of HWIL simulation has been described in reference [16]. The 6DoF model [19] runs in the HiGale real-time simulator, which is a multi-model and multi-target real-time simulation platform that runs on National Instruments' (NI) high-performance PXI (Peripheral Component Interconnection extensions for Instrumentation) controllers. HiGale is seamlessly connected with MATLAB, enabling Simulink to run in real time. The hardware interface is an NI multifunction reconfigurable I/O device, which is programmed for digital I/O, RS-422, and IMU pulses. The 6DoF model of a hypersonic vehicle includes 6DoF dynamics and kinematics models, aerodynamic models, mass/inertia models, Earth models, engine models, and guidance and control system models. The 6DoF model is the input source of the SINS in the HWIL simulation. The organic fusion of the 6DoF model and the SINS/GPS will enable the hypersonic vehicle guidance navigation and control system to conduct evaluation tests in the HWIL simulation system. On the basis of reference [16], a GPS simulator is added to the simulation platform to form a SINS/GPS integrated navigation system with the original inertial navigation system. Satellite-related data generated by GPS simulators can be used to verify the loose and tight coupling integrated SINS/GPS navigation simulation algorithms.
The HWIL simulation navigation simulator supports an IMU simulator and IMU model, a GPS simulator and GPS model, and a servo simulator (digital servo) and servo model, all in the loop mode. By using the IMU model and GPS model, system portability can be realized, which is highly beneficial for field test verification. These two models are described in detail below.

Inertial Measurement Unit (IMU) Model
The position, velocity, and attitude angle of the 6DoF model in the launch-centered Earth-fixed (LCEF) frame are the navigational theoretical parameters of a hypersonic vehicle. The LCEF frame is defined as: the coordinate origin is fixed to the launch point. The x-axis is in the horizontal plane of the launch point and points to the target direction; the y-axis is perpendicular to the horizontal plane of the launch point and points upward; and the z-axis and the x, y-axis form the right-handed coordinate system. The LCEF frame is usually represented by the g frame.
Details of the theoretical input of specific force and angular velocity, the IMU error model, and the velocity increment and angular increment of specific force and angular velocity in the 6DoF model are explained below. Figure 1 indicates the location of the IMU model in the HWIL simulation system. Figure 2 presents the implementation of the IMU model. The specific force and angular velocity are obtained from a high-precision 6DoF model. Based on these input parameters, an IMU error model is included. Finally, the continuous IMU signal is subjected to integral sampling and quantization to obtain the IMU pulses.

Theoretical Input of Specific Force and Angular Velocity
The specific force b f , which is the theoretical output of the accelerometer, is the sensitive force acting upon a unit mass other than gravity in the inertial frame; it is expressed as follows [16]: Equation (1) indicates that the specific force obtained from the 6DoF model is a result of the combined action of the engine thrust P , aerodynamic force R , actuator control force F c , and additional Coriolis force k ' F of flight vehicles. The Coriolis force is a description of the offset of the

Theoretical Input of Specific Force and Angular Velocity
The specific force b f , which is the theoretical output of the accelerometer, is the sensitive force acting upon a unit mass other than gravity in the inertial frame; it is expressed as follows [16]: Equation (1) indicates that the specific force obtained from the 6DoF model is a result of the combined action of the engine thrust P , aerodynamic force R , actuator control force F c , and additional Coriolis force k ' F of flight vehicles. The Coriolis force is a description of the offset of the

Theoretical Input of Specific Force and Angular Velocity
The specific force f b , which is the theoretical output of the accelerometer, is the sensitive force acting upon a unit mass other than gravity in the inertial frame; it is expressed as follows [16]: Equation (1) indicates that the specific force obtained from the 6DoF model is a result of the combined action of the engine thrust P, aerodynamic force R, actuator control force F c , and additional Sensors 2020, 20, 5418

of 24
Coriolis force F k of flight vehicles. The Coriolis force is a description of the offset of the linear motion of the mass point in the rotating system due to inertia relative to the linear motion of the rotating system, which comes from the inertia of the object motion [20]. The specific force reflects the actual motion state of a vehicle's centroid, which is different from that in the traditional trajectory generator.
In the 6DoF model, the angular velocity ω b in the launch-centered inertial (LCI) frame is the theoretical output of the gyroscope. Equation (2) suggests that the angular velocity generated by the 6DoF model is a result of the combined action of all types of moment and it reflects the actual angular motion of the vehicle revolving about its centroid; this is different from the traditional trajectory generators. In the LCI frame, the dynamic equation of the centroid of the hypersonic vehicle is expressed as [20].
where I is the inertia tensor of the vehicle, M st is the aerodynamic stabilizing moment, M c is the control moment, M d is the damping moment, M rel is the additional relative moment, and M k is the additional Coriolis moment. The above moments is defined as follows [21]: during the flight of the aircraft, the point of action of the aerodynamic force does not coincide with the center of mass, so the aerodynamic force will form a rotational moment on the center of mass, and this moment is called the aerodynamic stability moment; by changing the thrust direction of the aerodynamic engine, the force and moment that control the flight of the aircraft are generated, and this moment is defined as the control moment; when the aircraft rotates relative to the atmosphere, the atmosphere will have a damping effect on it, and the acting moment is defined as the damping moment; the relative force and Coriolis inertial force generated by the relative flow of fuel in the aircraft create moments on the center of mass, which are defined as additional relative moment and the additional Coriolis moment. The detailed calculation formula can refer to reference [21] 3.

Lever Arm Effects
The lever arm effect is a phenomenon, which means that when the rigid carrier has angular motion relative to the inertial space, the inertial sensors of the two inertial navigation systems at different positions will measure different specific forces, so that the calculated speed and position are also different [22]. In inertial navigation, the difference information of the main and sub inertial navigation system navigation parameters caused by the lever arm effect has nothing to do with the error propagation characteristics of the sub inertial navigation system. If it is not eliminated, it will affect the estimation accuracy of the error parameters of the sub-inertial navigation system and then affect the navigation performance [23].
The schematic diagram of the lever arm effect is shown in Figure 3, ib is the rotational angular velocity of the carrier relative to the inertial space. The accelerometer is installed at the fixed point P in the carrier frame; R O is the position vector of the origin of the carrier frame relative to the origin of the inertial frame; R P is the position vector of the point P relative to the origin of the inertial frame; r P is the position vector of the point P relative to the origin of the carrier frame.
Ideally, the installation point of the accelerometer should be located on the center of the carrier, that is the origin of the carrier frame O b , and when the mounting point P of the accelerometer deviates from the center of the carrier (r P is not zero). The lever arm effect error δ f caused by r P can be expressed as: where the subscript i represents the differential between the relative and inertial frame; the subscript b represents the differential between the relative and the carrier frame. is the specific force measured by the accelerometer at point P.
Sensors 2020, 20, x FOR PEER REVIEW 6 of 25 Figure 3. Schematic diagram of lever arm effect.
Ideally, the installation point of the accelerometer should be located on the center of the carrier, that is the origin of the carrier frame Ob, and when the mounting point P of the accelerometer deviates from the center of the carrier ( r P is not zero). The lever arm effect error f δ caused by r P can be expressed as: where the subscript i represents the differential between the relative and inertial frame; the subscript b represents the differential between the relative and the carrier frame.
represents the linear velocity of the origin of the carrier frame relative to the inertial frame; represents the linear velocity of the point P relative to the inertial frame; represents the linear velocity of the point P relative to the carrier frame.
is the specific force measured by the accelerometer at point Ob, and is the specific force measured by the accelerometer at point P.

Gyroscope and Accelerometer Error Models
The gyroscope error model is given as follows: is the gyroscope bias vector; b ω is the angular velocity vector of the gyroscope; g M is

Gyroscope and Accelerometer Error Models
The gyroscope error model is given as follows: where B b g is the gyroscope bias vector; ω b is the angular velocity vector of the gyroscope; M g is the gyroscope misalignment error; ε g is the gyroscope random noise vector; and S gi is the gyroscope scale factor error, i = x, y, z.
The accelerometer error model is given as follows: where B b a is the accelerometer bias vector; f b is the accelerometer specific force; M a is the accelerometer misalignment error; D a is the error associated with the acceleration quadratic term; ε a is the accelerometer random noise vector; and S ai is the accelerometer scale factor error, i = x, y, z.

Implementation of Integral Quantification
The specific force and angular velocity of the IMU are output in the form of velocity increment and angular increment, which are the integral results of the accelerometer theoretical value f b and the Sensors 2020, 20, 5418 7 of 24 gyroscope theoretical value ω b in unit time, respectively, as shown in Equations (8) and (9). The velocity increment and angular increment are further quantified and output as pulse numbers. Inside the IMU, the output quantification is generally performed after the high-frequency sampling process; that is, the IMU has the characteristics of internal high-frequency sampling and external low-frequency incremental output. In Figure 4, the acceleration output is taken as an example to illustrate the method of implementation of the velocity increment and the number of pulse.
accelerometer misalignment error; a D is the error associated with the acceleration quadratic term; a ε is the accelerometer random noise vector; and ai S is the accelerometer scale factor error, i = x, y, z .

Implementation of Integral Quantification
The specific force and angular velocity of the IMU are output in the form of velocity increment and angular increment, which are the integral results of the accelerometer theoretical value b f and the gyroscope theoretical value b ω in unit time, respectively, as shown in Equations (8) and (9). The velocity increment and angular increment are further quantified and output as pulse numbers. Inside the IMU, the output quantification is generally performed after the high-frequency sampling process; that is, the IMU has the characteristics of internal high-frequency sampling and external low-frequency incremental output. In Figure 4, the acceleration output is taken as an example to illustrate the method of implementation of the velocity increment and the number of pulse. As seen in Figure 4, through the sampling module, Simulink is divided into two parts, namely, 1 ms cycle period and 5 ms cycle period. For the 6DoF model, the 1 ms cycle is the simulation period (a smaller simulation period can be used according to the simulation accuracy need), and the 5 ms cycle is the running cycle of the simulation inertial device. Through the Simulink integral module, in the 1 ms cycle period, the acceleration/angular velocity information is accumulated every 1 ms within 5 ms, and the internal high-frequency sampling of the inertial device is simulated, so that there is no loss of the acceleration/angular velocity high-frequency information. In the 5 ms cycle period, the incremental output information of the inertial device is simulated by the Simulink delay block and the subtraction block. The incremental information of the IMU is quantified to obtain the quantified pulse number, and the implementation equations are obtained as follows: As seen in Figure 4, through the sampling module, Simulink is divided into two parts, namely, 1 ms cycle period and 5 ms cycle period. For the 6DoF model, the 1 ms cycle is the simulation period (a smaller simulation period can be used according to the simulation accuracy need), and the 5 ms cycle is the running cycle of the simulation inertial device. Through the Simulink integral module, in the 1 ms cycle period, the acceleration/angular velocity information is accumulated every 1 ms within 5 ms, and the internal high-frequency sampling of the inertial device is simulated, so that there is no loss of the acceleration/angular velocity high-frequency information. In the 5 ms cycle period, the incremental output information of the inertial device is simulated by the Simulink delay block and the subtraction block. The incremental information of the IMU is quantified to obtain the quantified pulse number, and the implementation equations are obtained as follows: where ∆V k is the integral increment in the 5 ms cycle period, corresponding to accIn in Figure 4; ∆V k−1 is the margin after the quantification of the last period pulse, with the margin magnitude being less than 1 pulse equivalent, corresponding to accIn_r in Figure 4; ∆V k is the margin after the quantization of the current beat pulse, corresponding to accOut_r in Figure 4; ∆V Pul is the pulse number after current quantification, corresponding to accPul in Figure 4; PULSE_ACC is the pulse equivalent.

Satellite Receiver Simulation of Integrated Navigation
In the HWIL simulation diagram (Figure 1), the location of the GPS receiver model in the HWIL simulation is indicated. Figure 5 presents the implementation of the GPS receiver model. We specify the GPS satellite constellation through a GPS broadcast ephemeris file. The daily GPS broadcast ephemeris file (brdc) is obtained by merging individual site navigation files. The GPS broadcast ephemeris files are then used to generate the simulated pseudo-range and Doppler shift for the visible GPS satellites. Through the parameters provided by the ephemeris, the corresponding satellite parameters such as satellite position, speed, and clock correction can be calculated.
In the HWIL simulation diagram (Figure 1), the location of the GPS receiver model in the HWIL simulation is indicated. Figure 5 presents the implementation of the GPS receiver model. We specify the GPS satellite constellation through a GPS broadcast ephemeris file. The daily GPS broadcast ephemeris file (brdc) is obtained by merging individual site navigation files. The GPS broadcast ephemeris files are then used to generate the simulated pseudo-range and Doppler shift for the visible GPS satellites. Through the parameters provided by the ephemeris, the corresponding satellite parameters such as satellite position, speed, and clock correction can be calculated.

Satellite Receiver Simulation of Integrated Navigation
As shown in Figure 6, the steps of the simulation of satellite data in the loose and tight coupling integrated navigation are as follows. First, the GPS satellite ephemeris is used to calculate the position and velocity of the satellite; then, the theoretical position and velocity of the hypersonic vehicle 6DoF model and the position and velocity of the satellite are used to calculate the satellite's pseudo-range and pseudo-range rate or Doppler shift, which are the simulation data required for tight coupling. Finally, using the combination of the position and velocity of the satellite, the position and velocity of the receiver are calculated by the least squares method, which are the simulation data required for loose coupling.

Satellite Receiver Simulation of Integrated Navigation
As shown in Figure 6, the steps of the simulation of satellite data in the loose and tight coupling integrated navigation are as follows. First, the GPS satellite ephemeris is used to calculate the position and velocity of the satellite; then, the theoretical position and velocity of the hypersonic vehicle 6DoF model and the position and velocity of the satellite are used to calculate the satellite's pseudo-range and pseudo-range rate or Doppler shift, which are the simulation data required for tight coupling. Finally, using the combination of the position and velocity of the satellite, the position and velocity of the receiver are calculated by the least squares method, which are the simulation data required for loose coupling.  Tight coupling SINS/GPS integrated navigation requires pseudo-range and pseudo-range rate measurements of the satellite. In the simulation platform, the pseudo-range and pseudo-range rate can be calculated from the position and velocity of the satellite and receiver. The calculation equation for the position and velocity of the satellite has been described in references [24,25].
The steps for the calculation of the pseudo-range and pseudo-range rate based on the position and velocity of the satellite and receiver are as follows: Here, ∆p =[∆x, ∆y, ∆z] is the position difference between the satellite and the receiver, [x u , y u , z u ] T is the position of the receiver in the Earth-centered Earth-fixed frame (ECEF, e frame), and [x sv , y sv , z sv ] T is the position of the satellite in the ECEF frame. The ECEF is defined as: the origin is the center of the Earth, and the x-and y-axes are in the equatorial plane of the Earth. The x-axis points to the prime meridian, while the z-axis is the Earth's rotation axis [26].
In Equation (12), ∆t tr is the transmission time of the satellite signal from the satellite to the receiver, c is the speed of light.
x sv , y sv , z sv T is the position coordinates of the satellite after considering the rotation effect of the earth, ω ie is the earth's rotation angular velocity.
∆p =[∆x , ∆y , ∆z ] T is the position difference between the satellite and the receiver after considering the rotation effect. Then, the approximate distance from the satellite to the receiver is: After considering the clock correction of the GPS, ionospheric error, and tropospheric error, the pseudo-range is obtained as: Here, T iono is the time delay caused by ionospheric errors, δt s is the satellite clock correction. The pseudo-range rate is obtained as: In Equation (17), v x , v y , v z is the velocity difference between the satellite and the receiver.

Loose Coupling Data Simulation
Loose coupling requires the position and velocity of the satellite receiver, which are obtained from the position and velocity parameters of the satellite and the pseudo-range and pseudo-range rate (or Doppler shift) measurements.
The position and velocity of the receiver are calculated by the least squares method. The specific calculation steps are as follows.
The mathematical model of single-mode pseudo-range measurements is: Equations (17) and (18) represents m pseudo-range measurements, [x u , y u , z u ] is the position of the receiver, b is the deviation between the local clock of the receiver and GPST (global positioning system time); and n ρ is the error in the pseudo-range measurement. The known quantities ρ i are pseudo-range measurements and satellite coordinates x s i , y s i , z s i .
The pseudo-range measurement can be linearized and written in the form of a matrix as follows: In Equation (19), δρ = [δρ 1 , δρ 2 , · · · , δρ m ] T , H = [u 1 , u 2 , · · · , u m ] T , n ρ = n ρ 1 , n ρ 2 , · · · , n ρ m T , u i = Equation (20) obtains the correction value between the initialization and the real point after the first power linearization; this correction value can be applied to update the initial point to obtain the corrected solution; that is: Then, x 1 can be used as the starting point to repeat the process from linearization to Equation (21) in order to obtain a new correction value dx 1 and the previous solution is updated.

Geometric Dilution of Precision
When performing satellite navigation solutions, the geometric dilution of precision can reflect the quality of satellite positioning to a certain extent. Therefore, the geometric dilution of precision is selected in the simulation platform as the index for evaluating the satellite signal generated by the satellite simulator [27,28].
According to the covariance matrix of the state error obtained by the least squares method in the previous section, the covariance matrix of the position error at this time can be obtained as follows: Here, δx u is used to distinguish from the updated dx i of each iteration in the previous section; R is the covariance matrix of the noise vector in the pseudo-range measurement. It is generally assumed that the measurement noises of different satellites are independent of each other, so R is a diagonal matrix, expressed as diag σ 2 1 , σ 2 2 , · · · , σ 2 m . Obviously, σ 2 i is the noise power, which measures the quality of the i-th satellite pseudo-range measurement? R = σ 2 I is obtained and I can be a unit matrix m × m. Substituting R = σ 2 I in Equation (22), we obtain: Let H T H −1 be written with h i.j , where h i.j represents the i-th row and j-column elements of the matrix. Hence, It can be seen from Equation (24)  The dilution of precision can be seen as a linear mapping from the measurement error to the state estimation error. In the case where the measurement errors are the same, a larger dilution of precision may result in a larger state estimation error, and a smaller dilution of precision results in a considerably smaller state estimation error. It can be seen from the definition of the dilution of precision that it is independent of the actual measurement noise, while it can only be related to H T H −1 , which is directly calculated from the H matrix.

Global Positioning System (GPS) Error
The flying height of hypersonic vehicles is more than 20 km. The GPS errors to be considered mainly include ionospheric errors and satellite clock corrections.

Ionospheric Error
The ionosphere is a layer of the Earth's atmosphere that is ionized by the solar ray; it is 50-1000 km from the surface of the Earth and the inner boundary of the Earth's magnetosphere [27][28][29]. Under the intense illumination of sunlight, the neutral gas molecules in the ionosphere are ionized, resulting in the generation of a large number of positive ions and free electrons. These positive ions and free electrons affect the propagation of navigational waves.
When the GPS satellite signal passes through the ionosphere, a certain delay exists, which is proportional to the number of free electrons encountered [30,31]. The electron density is affected by factors such as local time, geomagnetic latitude, and sunspot cycle. The ionospheric correction model used is the Klobuchar ionospheric delay model [32]. The Klobuchar model abstracts ionospheric delay as a cosine function. A is the amplitude of the cosine curve; PER is the period of the cosine curve. The composition of the ionospheric correction algorithm for the GPS receiver requires the receiver's approximate latitude and longitude positions (φ U , λ U ) as well as the azimuth and elevation angles (Az, El) relative to each satellite.
The calculation steps are shown in Figure 7.
The ionospheric correction model used is the Klobuchar ionospheric delay model [32]. The Klobuchar model abstracts ionospheric delay as a cosine function. A is the amplitude of the cosine curve; PER is the period of the cosine curve. The composition of the ionospheric correction algorithm for the GPS receiver requires the receiver's approximate latitude and longitude positions as well as the azimuth and elevation angles ( ) Az, El relative to each satellite.
The calculation steps are shown in Figure 7.   In Figure 7, ψ is the center angle of the Earth, φ I , λ I , and φ m are the geodetic latitude, geodetic longitude, and geomagnetic latitude of the earth projection of the ionospheric intersection point, respectively, t is the local time of the ionospheric intersection point, F is the ionospheric mapping function, T iono is the ionospheric delay. When A is smaller than 0, A equal to 0 is taken, when PER is smaller than 72,000, PER equal to 72,000 is taken, α i , β i is the ionospheric parameters of the i-th satellite in the satellite navigation message.

Satellite Clock Correction
A precise clock is critical in GPS. The GPS receiver receives signals from satellites whose transmission time is proportional to the pseudo-range measurement. Then, the relationship between the transmission time and pseudo-range measurement is: Here, t r is the propagation time of the signal from the satellite to the receiver, ρ is the pseudo-range measurement, and c is the speed of light.
The satellite emission signal time can be expressed as: Here, t sv is the satellite signal transmission time and t m is the signal reception time.
The satellite clock error can be calculated from the polynomial factors obtained from the GPS navigation message.
∆t sv = a f 0 + a f 1 (t−t oc ) + a f 2 (t−t oc ) 2 Here, ∆t sv is the satellite clock correction value, a f 0 is the satellite clock error, a f 1 is the small error, a f 2 is the small frequency drift, t oc is the clock reference point, and the value of (t−t oc ) should be corrected within seconds of a week.
Then, the correction value of the signal propagation time is: Here, ∆t r is a relativistic correction value, which can be expressed as: Here, e is the eccentricity of the satellite orbit, a is the semi-major axis length, and E k is the eccentric anomaly.
The satellite clock difference can be defined as: The correction equation for the pseudo-range measurement is: In the above parameters,a f 0 ,a f 1 ,a f 2 , e, a, and E k are obtained from the broadcast ephemeris.

Simulation Verification
The simulation verification used the 1100 s trajectory data of the hypersonic vehicle to verify the simulation. The initial state of the trajectory is shown in Table 1. The typical trajectory of a hypersonic boost-glider is presented in Figure 8 [18,33].  The typical trajectory of a hypersonic boost-glider is presented in Figure 8 [18,33].

IMU Data Simulation
The relationship between the number of accelerometer pulses and the specific force is depicted in Figure 9.
is the output of the accelerometer pulse number. Due to the quantification process, the number of pulses is an integer; is the accelerometer specific force. Accelerometer

IMU Data Simulation
The relationship between the number of accelerometer pulses and the specific force is depicted in Figure 9.

IMU Data Simulation
The relationship between the number of accelerometer pulses and the specific f in Figure 9.
is the output of the accelerometer pulse number. Due to th process, the number of pulses is an integer; is the accelerometer specific force pulse equivalent is 0.0001 2 m / s . The relationship between the number of gyroscope angular velocity is depicted in Figure 10.
is the output of the gyroscope pu is the gyroscope angular velocity. Gyroscope pulse equivalent × 0.1 (π / 180) / is the output of the accelerometer pulse number. Due to the quantification process, the number of pulses is an integer;

IMU Data Simulation
The relationship between the number of accelerometer pulses and the specific force is depicted in Figure 9.
is the output of the accelerometer pulse number. Due to the quantification process, the number of pulses is an integer; is the accelerometer specific force. Accelerometer pulse equivalent is 0.0001 2 m / s . The relationship between the number of gyroscope pulses and the angular velocity is depicted in Figure 10.
is the output of the gyroscope pulse number and is the gyroscope angular velocity. Gyroscope pulse equivalent × 0.1 (π / 180) / 3600 rad / s . is the accelerometer specific force. Accelerometer pulse equivalent is 0.0001 m/s 2 . The relationship between the number of gyroscope pulses and the angular velocity is depicted in Figure 10.

yaw angle 0°
The typical trajectory of a hypersonic boost-glider is presented in Figure 8

IMU Data Simulation
The relationship between the number of accelerometer pulses and the specific f in Figure 9.
is the output of the accelerometer pulse number. Due to th process, the number of pulses is an integer; is the accelerometer specific force pulse equivalent is 0.0001 2 m / s . The relationship between the number of gyroscope angular velocity is depicted in Figure 10.
is the output of the gyroscope pu is the gyroscope angular velocity. Gyroscope pulse equivalent × 0.1 (π / 180) / is the output of the gyroscope pulse number and pitch angle 90° roll angle 0° yaw angle 0° The typical trajectory of a hypersonic boost-glider is presented in Figure 8 [18,33].

IMU Data Simulation
The relationship between the number of accelerometer pulses and the specific force is depicted in Figure 9.
is the output of the accelerometer pulse number. Due to the quantification process, the number of pulses is an integer; is the accelerometer specific force. Accelerometer pulse equivalent is 0.0001 2 m / s . The relationship between the number of gyroscope pulses and the angular velocity is depicted in Figure 10.

yaw angle 0°
The typical trajectory of a hypersonic boost-glider is presented in Figure 8 [18,33].

IMU Data Simulation
The relationship between the number of accelerometer pulses and the specific force is depicted in Figure 9.
is the output of the accelerometer pulse number. Due to the quantification process, the number of pulses is an integer; is the accelerometer specific force. Accelerometer pulse equivalent is 0.0001 2 m / s . The relationship between the number of gyroscope pulses and the angular velocity is depicted in Figure 10.
is the output of the gyroscope pulse number and is the gyroscope angular velocity. Gyroscope pulse equivalent × 0.1 (π / 180) / 3600 rad / s .

GPS Receiver Data Simulation
The data of the trajectory generator are imported into the GPS simulator to obtain corresponding GPS satellite data. Since the attitude of the hypersonic vehicle changes constantly during the actual flight, it may affect the satellite receiver's ability to determine the visible satellites. As shown in Figure 11, the number of visible satellites varies.

GPS Receiver Data Simulation
The data of the trajectory generator are imported into the GPS simulator to obtain corresponding GPS satellite data. Since the attitude of the hypersonic vehicle changes constantly during the actual flight, it may affect the satellite receiver's ability to determine the visible satellites. As shown in Figure 11, the number of visible satellites varies. Figure 10. Output of gyroscope pulse output and gyroscope angular velocity.

GPS Receiver Data Simulation
The data of the trajectory generator are imported into the GPS simulator to obtain corresponding GPS satellite data. Since the attitude of the hypersonic vehicle changes constantly during the actual flight, it may affect the satellite receiver's ability to determine the visible satellites. As shown in Figure 11, the number of visible satellites varies. The positional dilution of precision of satellites is shown in Figure 12. The positional dilution of precision of satellites is shown in Figure 12. Satellite data can be divided into tight and loose coupling data. Tight coupling data include satellite pseudo-range and pseudo-range rate (Doppler shift); Figure 13 shows the variation of these parameters. Satellite data can be divided into tight and loose coupling data. Tight coupling data include satellite pseudo-range and pseudo-range rate (Doppler shift); Figure 13 shows the variation of these parameters.
The data output from the satellite simulator is input to the satellite receiver for verification, and the position and velocity errors between the original trajectory and the calculated trajectory are obtained as shown in Figures 14 and 15.
It can be seen from the simulation results that the position errors are within 20 m (horizontal direction) and velocity errors are within 0.  Satellite data can be divided into tight and loose coupling data. Tight coupling data include satellite pseudo-range and pseudo-range rate (Doppler shift); Figure 13 shows the variation of these parameters. satellite pseudo-range and pseudo-range rate (Doppler shift); Figure 13 shows the variation of these parameters.
(a) (b) Figure 13. The variation of tight coupling data: (a) satellite pseudo-range variation; (b) satellite Doppler shift.
The data output from the satellite simulator is input to the satellite receiver for verification, and the position and velocity errors between the original trajectory and the calculated trajectory are obtained as shown in Figures 14 and 15.

Strapdown Inertial Navigation Systems (SINS)/GPS Receiver Data Simulation
In the simulation verification of the integrated navigation algorithm, the gyroscope constant drift [34] is 3/h(1σ), the gyroscope random error is 0.3/ √ h(1σ), the gyroscope scale factor error is 100 × 10 −6 , the accelerometer constant bias is 1 × 10 −3 g(1σ); the accelerometer random error is 1 × 10 −4 g/ √ Hz(1σ), the accelerometer scale factor error is 100 × 10 −6 ; the initial roll angle, yaw, pitch errors are 60 , 20 , 20 , the initial velocity error is 0.05 m/s, the initial position is 5 m [18,33]. Figures 16-30 present the simulation results of the SINS/GPS integrated navigation system conducted 50 times. The difference between them is that random errors are added to the gyroscope constant drift and accelerometer constant bias to better simulate the real inertial navigation system. Note that the three attitude angle errors of the integrated navigation system can basically converge to within 0.05 • ; the speed error in three directions can reach 0.2 m/s; and the position error of the direction is within 10 m; the gyroscope constant drift can converge to 3/h; the accelerometer constant bias can converge to 1 × 10 −3 g. This result meets the navigation requirements of hypersonic vehicles.

Time Analysis
The results of the HiGale real-time simulator real-time simulation are shown in the Figure 31. In the figure, TotalRunTimer is the system running time, the unit is second, StepTimer is the HiGale simulation step size, the unit is microsecond, and Turnaround is the simulation software time consumption, the unit being the microsecond. This article uses a 1 ms simulation cycle, but HiGale internally uses the fourth Runge-Kutta method, that is, using a 500 µs cycle, Turnaround is the simulation software running time, both are less than 500 µs. The results of the HiGale real-time simulator real-time simulation are shown in the Figure 31. In the figure, TotalRunTimer is the system running time, the unit is second, StepTimer is the HiGale simulation step size, the unit is microsecond, and Turnaround is the simulation software time consumption, the unit being the microsecond. This article uses a 1 ms simulation cycle, but HiGale internally uses the fourth Runge-Kutta method, that is, using a 500 μs cycle, Turnaround is the simulation software running time, both are less than 500 μs .

Conclusions
The simulation platform for the inertial navigation algorithm of hypersonic vehicles based on flight mechanics proposed in this paper is different from the traditional trajectory generator. It can generate a real aircraft trajectory, reproduce the real flying environment of the aircraft in the sky as realistically as possible, and can be coordinated with the flight control system in the HWIL simulation to verify and evaluate the performance index of the flight control system. The simulation platform can also generate simulation data to verify the accuracy of the integrated navigation algorithm. The simulation test in this paper verifies the correctness of the design method.
Author Contributions: Conceptualization, K.C. and F.S.; methodology, K.C. and F.S.; software, F.S.; validation, K.C. and J.Z.; data curation, F.S.; writing-original draft preparation, K.C. and X.W.; writing-review and editing, K.C. and X.W. All authors have read and agreed to the published version of the manuscript.

Conclusions
The simulation platform for the inertial navigation algorithm of hypersonic vehicles based on flight mechanics proposed in this paper is different from the traditional trajectory generator. It can generate a real aircraft trajectory, reproduce the real flying environment of the aircraft in the sky as realistically as possible, and can be coordinated with the flight control system in the HWIL simulation to verify and evaluate the performance index of the flight control system. The simulation platform can also generate simulation data to verify the accuracy of the integrated navigation algorithm. The simulation test in this paper verifies the correctness of the design method.