SINS/BDS Integrated Navigation for Hypersonic Boost-Glide Vehicles in the Launch-Centered Inertial Frame

According to the trajectory specialty of hypersonic boost-glide vehicles, a strapdown inertial navigation system/BeiDou navigation satellite system (SINS/BDS) algorithm based on the launch-centered inertial (LCI) frame for hypersonic vehicles is proposed. First, the related frame system, especially the launch earth-centered inertial (LECI) frame, and the SINS mechanization in the LCI frame are introduced. Second, SINS discrete updating algorithms in the LCI frame for the compensation of coning, sculling, and scrolling effects are deduced in the attitude, velocity, and position updating algorithms, respectively. Subsequently, the Kalman filter of the SINS/BDS integrated navigation in the LCI frame is obtained. )e method of converting BDS receiver position and velocity from the Earth-centered Earth-fixed (ECEF) frame to the LCI frame is deduced through the LECI frame. Finally, taking the typical hypersonic boost-glide vehicles as the object, the SINS/BDS algorithm vehicle field test and hardware-in-the-loop simulation are performed.


Introduction
A hypersonic vehicle is an aircraft that has the ability to fly and cruise for a long time at more than five times the speed of sound in near space; it is also called a "near-space hypersonic vehicle" [1]. e near space is the atmospheric space 20-100 km from the ground, and it is located below the low-orbit satellite operation area and above the traditional aircraft flight area. It has great military and civilian application value [2,3].
Scholars from all over the world have studied many different navigation frames for different hypersonic vehicle navigation algorithms [4]. Stephen studied the strapdown algorithm in the earth-centered inertial frame for the SHEFEX2 hypersonic vehicle [5]. Mu et al. studied an integrated navigation algorithm of hypersonic aircraft in the ECI frame [6]. Wei et al. studied a hypersonic vehicle GFSINS/GPS/CNS algorithm using the ECEF frame [7]. Yu et al. used the NEU geography coordinate system as the hypersonic vehicle navigation reference system [8]. e X-43A used an INS/GPS integrated navigation system and the geography coordinate system navigation algorithm [9]. Kang et al. studied an INS/BDS integrated navigation method in the ENU geography coordinate system [10]. Li studied the aerial alignment method of the integrated navigation system in the local horizontal coordinate system [11].
Because a hypersonic boost-glide vehicle is launched vertically, its trajectory is similar to that of the launch vehicle. If the local horizontal frame is used as the navigation coordinate system of the hypersonic boost-glide vehicle, the attitude angle will be singular during the navigation calculation process. erefore, the LCI frame is often used as the navigation coordinate system of the launch vehicle, which will not cause the attitude angle to be singular, making the LCI frame more advantageous over the local horizontal frame [4,12,13]. According to the navigation solutions in the LCI frame, scholars from all over the world have designed some hypersonic vehicle flight control algorithms. Zhang et al. studied stellar-inertial integrated guidance to responsive launch vehicles in an LCI frame [12]. Zhao et al. built an INS/CNS/GPS system in the LCI frame for vehicles for reducing the effects of an imprecise Earth model and completing long-duration maneuvers within a wide flight envelope [13]. Our group proposed an algorithm of SINS in the LCI frame for hypersonic boost-glide vehicles and presented a solution in the local horizontal frame from the LCI frame.
In the near-space hypersonic flight, factors such as aerodynamic disturbance, scramjet operation, and steering gear operation affect the accuracy of the SINS of hypersonic vehicles. e Chinese BDS navigation system has become a world-class global navigation satellite system (GNSS) [14,15]. e application of BDS is less in hypersonic vehicles. In the present study, based on the high-dynamic characteristics of hypersonic vehicles, we propose updating algorithms for the attitude, velocity, and position of strapdown inertial navigation as well as a SINS/BDS integrated navigation system in the LCI frame. In addition, we conduct a test on the SINS/BDS integrated navigation vehicle and hardware-in-the-loop (HWIL) simulation of the LCI frame.

SINS Mechanization in Launch-Centered
Inertial Frame

Common Frame Systems.
In the SINS, the commonly used frame systems include the Earth-centered inertial frame (ECI, i frame), Earth-centered Earth-fixed frame (ECEF, e frame), launch Earth-centered inertial frame (LECI, t frame), launch-centered Earth-fixed frame (LCEF, g frame), launchcentered inertial frame (LCI, a frame), and body-fixed frame (BF, b frame) [4]. is study uses the frame systems described below and the LCI frame as the navigation frame system.

Relationship between
where φ, ψ, and c are the vehicle pitch, yaw, and roll angles in the LCI frame, respectively. R i (θ) is the DCM rotating around the i-axis (i � x, y, z) by angle θ.

LCEF Frame and ECEF Frame.
e DCM from the ECEF frame to the LECF frame is R g e , as shown in where ϕ 0 , λ 0 , and α 0 are the initial latitude, initial longitude, and launch azimuth of the vehicle. In addition, at the time of launch, the coordinates of the LCEF, and LCI frames are in the same direction, and the coordinates of the LECI and ECEF frames are in the same direction. erefore, the DCMs R a t from the LECI frame to the LCI frame and R g e are the same. During flight, R a t and R g e are constant values.

Initial Values of the Launch Point.
e launch point position in the ECEF frame is P e 0 � [x e0 , y e0 , z e0 ] T , as given by where H 0 is the initial height, N 0 � a e / ���������� 1 − e 2 sin 2 L 0 , e is the eccentricity, and a e is the semimajor axis.
At the time of the launch, the coordinates of the LECI and ECEF frames point in the same direction, so P t 0 � P e 0 .

Geocentric Vector of the Launch Point.
e geocentric vector of the launch point in the LCEF frame is P g 0e : At the time of launch, the coordinates of the LCEF and LCI frames point in the same direction, so P a 0e � P g 0e .

Rotational Angular Velocity in the LCI Frame.
e Earth's rotation angular velocity in the LECF frame relative to the LCI frame is ω g ag .

SINS Mechanization in the LCI Frame.
e navigation equation of the SINS in the LCI frame is given by the following equation, and its mechanization is shown in Figure 1: As shown in (6) and Figure 1, the SINS algorithm is executed by the navigation computer, and steps involved are as follows. First, the vehicle attitude is obtained by integrating angular velocity ω b measured by gyroscopes, and then attitude matrix R a b is obtained. Second, specific force f b measured by accelerometers is converted to the LCI frame through attitude matrix R a b ; considering the effect of gravity g a , velocity V a of the vehicle is obtained by integrating the sum of f a and g a in the LCI frame. Finally, the position P a of the vehicle is obtained by integrating velocity V a , and P g � [x, y, z] T . ese three integral processes are termed SINS attitude updating, velocity updating, and position updating, respectively. e formula for calculating gravity in the LCI frame is found in reference: 2 Mathematical Problems in Engineering where g r ′ is calculated by (8) and g ω e is calculated by (9). r 0 is the unit vector of r, and r is calculated by (10): r � P a 0e + P a .
In these equations, μ is the coefficient of gravity; J � (3/2)J 2 , where J 2 is the harmonic coefficient; r � ���������� x 2 + y 2 + z 2 ; ϕ is the Earth-centered latitude, ϕ � arcsin(rw/(|r| * ω ie )); and rw � r · ω e . e initial value of position P a is P a 0 � [0, 0, 0] T , and the initial value of attitude matrix R b a is determined by the initial attitude angles φ 0 , ψ 0 , and c 0 in the initial alignment. e initial value V a 0 is obtained as follows:

SINS Numerical Updating Algorithm in the LCI Frame
e hypersonic flight is generally affected by coning, sculling, and scrolling effects [16,17]. erefore, the compensation algorithms for the three effects are derived in attitude updating [18], velocity updating [19], and position updating [19], respectively.
In this study, the two-sample updating algorithm is used; in other words, the angular velocity vector ω b from the gyroscope and the specific force measurement vector f b from the accelerometer are both assumed to be linear models; that is, where a, b, A, B are constant vectors; then, we can get the angular and velocity increment expressions, If the IMU is sampled twice at equal intervals in [t k−1 , t k ], the sampling time is t k−1/2 and t k , respectively, T � t k − t k−1 and h � T/2 are recorded; then we can get where Δθ 1 and ΔV 1 are the angular and velocity increments in the time period [t k−1 , t k−1 + h], and Δθ 2 and ΔV 2 are the angular and velocity increments in the time period

Attitude Updating Algorithm.
e quaternion attitude updating algorithm based on rotation vector is completed in two steps: (1) calculation of a rotation vector, which describes the state change of the vehicle attitude, and (2) update of quaternions. e recurrence formula of attitude quaternion updating is as follows: where the time before and after attitude updating is t k−1 and t k , during which the quaternions of attitude are q a b(k−1) and is calculated by angular velocity ω b ab (t), which is angular velocity ω b . Note that Φ k is the equivalent rotation vector of BF frame relative to the LCI frame, and the corresponding equivalent vector differential equation approximated as follows: Since the output of the optical gyroscope is an angle increment, it is necessary to solve the rotation vector by the Taylor series expansion method. Take Φ k as the equivalent rotation vector in the time period [t k−1 , t k ]; then, the Taylor Navigation computer . b Figure 1: Mechanization of SINS in the LCI frame.

Mathematical Problems in Engineering
From (12), the derivative of each order of the vehicle angular velocity is obtained as follows: From (13), the derivative of each order of the angular increment is obtained as follows: Further, because the attitude updating period T is generally in milliseconds, Φ k can also be regarded as a small magnitude.
erefore, when calculating the derivative of Φ k (τ) at τ � 0, according to (16), replace the second term of Φ k (τ) with the angle increment; then, we can get us, we can get the following equation from (16): Considering (18) and (19), the derivative of each order of (21) is obtained: According to (18) and (19), τ � 0 is substituted into (22); then, we have Substituting (23) into (19), we can get where Φ k (0) is the rotation vector in the time period T.
Since the time interval is 0, Φ k (0) � 0; therefore, From (14), the model coefficients a and b of the linear fitting of the vehicle angular velocity can be obtained by inverse solution; that is, 4 Mathematical Problems in Engineering Substituting (26) into (24), we can get To solve the rotation vector according to (27), the angle increments Δθ 1 and Δθ 2 in the two time periods are used. erefore, this formula is called the two-sample updating algorithm of the rotation vector.
erefore, the solution of (16) by using the two-sample algorithm is as follows: Rotation vector Φ k and its corresponding quaternion are related as follows: is substituted into (16) to complete the attitude updating.

Velocity Updating Algorithm.
When a carrier operates in the environment of line vibration and angular vibration, the sculling error may occur, which should not be ignored in a high-precision SINS algorithm. e specific force given by (6) indicates the time parameters of each component as follows: From (31), where V a k and V a k−1 are velocities of the LCI frame at t k and t k−1 , respectively, and are recorded as where ΔV a sf(k) and ΔV a g(k) are the velocity increments of specific force and gravity in time period T, respectively. e following recurrence forms can be obtained: e numerical integration algorithms for ΔV a sf(k) and ΔV a g(k) are discussed below: (1) First, we have the gravity velocity increment ΔV a g(k) . Even for hypersonic vehicles, the change in the gravitational vector in a short time period T is very small, so the integrand function of ΔV a g(k) is generally considered to be a slow time variable, which can be replaced with the value at t k−1/2 � (t k−1 + t k )/2, and ΔV a g(k) is approximated as follows: Since the navigation velocity and position at time t k are unknown, g a k−1/2 needs to be calculated by linear extrapolation: (2) Second, we have the specific force velocity increment ΔV a sf(k) . e integrand matrix of (33) is transformed to the matrix chain multiplication decomposition as follows. Hence, For t k−1 ≤ t ≤ t k , the relationship between the coordinate transformation matrix and the equivalent rotation vector is as follows: Equation (42) can be approximated as (Φ×)(Φ×) can be regarded as a second-order small quantity. Equations (20) and (40) can be approximated as By substituting (42) into (39), we obtain Mathematical Problems in Engineering Hence, From (14), the following formula can be obtained: Equations (44) and (47) are substituted into (43), when t � t k+1 , ΔV(t k−1 ) � 0, and Δθ(t k−1 ) � 0 and when t � t k , ΔV(t k ) � ΔV k , and Δθ(t k ) � Δθ k . Hence, where Here, ΔV rotk is the compensation quantity of the velocity rotation effect, caused by the rotation of the linear motion direction of the carrier.
ΔV sculk is the compensation term of the velocity sculling effect when the carrier has linear and angular vibrations. Equations (12) and (13) are substituted into (50) and integrated to obtain

Mathematical Problems in Engineering
From (14), By substituting (26) and (52) into (51), the two-sample algorithm for the compensation of the velocity sculling error is obtained as follows: us, velocity updating algorithm derivation is completed.

Position Updating Algorithm.
When the hypersonic vehicle is in sculling flight, a scrolling effect error is present in the solution of the SINS position. According to (35) and (36), where

By integrating (54) in time period
From (38), we obtain Using (48), the following can be calculated: Note that c 1 � 1/2 e above integral is obtained by a step integral method: us, Mathematical Problems in Engineering e above equation is substituted into (58) to obtain at is, the second integral increment of specific force is obtained.
Equation (63) is the compensation of the rotation effect in position calculation.
Equation (64) is the compensation of the scrolling effect in position calculation, and it indicates that the factors affecting the scrolling effect are the sculling effect and the coupling effect between the angular and linear motions of the vehicle.

S b
Δvk � T us, position updating algorithm derivation is completed.

Attitude Error Equation.
For the SINS, the transformation matrix error from the BF frame to the calculated frame is caused by the angular rate error between the two coordinate systems when they are rotated. e attitude mechanization in the LCI frame is as follows.
Considering the measurement and calculation errors, the transformation rate of the calculated transformation matrix is given as follows.
ere is an error angle φ a � ϕ x ϕ y ϕ z T between the actual navigation frame and the calculated frame, and Ψ a is the antisymmetric matrix of φ a :

Mathematical Problems in Engineering
Derivation of both sides of (74) gives the following: (76) By contrast, the differential equation of (70) is obtained as Comparing (76) and (77), we get Equation (78) is written in a vector form as e gyroscope measurement error is δω b ab � δω b ai + δω b ib , and because δω b ai � 0, δω b ib is assumed to be composed of the random constant drift ε b of the equivalent gyroscope. us, erefore, the attitude error equation is written in the component form as

Velocity Error Equation.
From (6), the linear equation of the first-order small disturbance for velocity can be obtained as e measurement error δf b of the accelerometer consists of the random constant zero bias ∇ b of the equivalent accelerometer. erefore, we can get the velocity error equation as follows.

Position Error Equation.
e position error equation under the LCI frame is

SINS/BDS Integrated Navigation Equation.
e state and measurement equation for SINS/BDS integrated navigation are given as follows.
where Z vp is shown in It is assumed that the measurement error δω b ab of the gyroscope is composed of the gyroscope's zero bias ε b and that the measurement error δf b of the accelerometer consists of the accelerometer's zero bias ∇ b . e state variable of SINS/BDS is adopted; that is, e state equation is Mathematical Problems in Engineering e 15-dimensional measurement matrix H vp is expressed as follows:

BDS Receiver Position and Velocity Conversion.
Since the BDS receiver outputs the position and velocity in the ECEF frame [20], they need to be converted to the LCI frame. Let the position and velocity of the BDS receiver in the ECEF frame be P e s and V e s , respectively. First, by converting these values in to the LECI frame, position P t s and velocity V t s in the LECI are obtained as shown in en, position P a s and velocity V a s of the BDS receiver in the LCI frame are obtained by converting position P t s and velocity V t s from the LECI frame to the LCI frame, as shown in

SINS/BDS Integrated Navigation Vehicle Field Test.
To test the feasibility of the SINS/BDS integrated navigation algorithm, a vehicle field test was conducted, as shown in Figure 2. e test was conducted at a high speed around Xi'an, and the track of the vehicle is shown in Figure 3. e initial latitude and longitude of the vehicle test were 34.19785°and 108.8284°, respectively, the initial height was 365.1 m, and the initial velocity was 0 m/s. e vehicle was equipped with a set of high-precision integrated navigation system as the reference, called the master INS. Its main performance indicators are as follows: position accuracy is 1.0 n mile/h, speed accuracy is 0.8 m/s, and attitude accuracy is 0.05°. e SINS/BDS used in this paper is called the slave INS. e gyroscope drift was 3°/h (1σ), the scale factor was 100 ppm, the variance of white noise was 0.3°/h (1σ), and the output period of gyroscope was 2 ms. e bias of the accelerometer was 1 × 10 −3 ·g (1σ), the scale factor was 100 ppm, the white noise variance was 1 × 10 −4 ·g (1σ), and the output period was 2 ms. e positioning accuracy of the BDS satellite receiver was 10 m (1 σ), the speed accuracy was 0.3 m/s (1σ), and the output period was 1 Hz. e master INS output the navigation solutions were in the local horizontal frame. e slave INS solutions were converted from the LCI frame to the local horizontal frame [21]. Both are then compared and made the difference as shown in  In the vehicle test, the azimuth angle error of SINS/BDS integrated navigation was approximately 0.05°, the pitch angle error was within 0.05°, the roll angle error was within 0.05°, the velocity error of east and north was within 0.2 m/s, and the position error was within 10 m.

SINS/BDS Integrated Navigation Test by HWIL
Simulation. Due to the limitation of vehicle velocity, the performance of SINS/BDS integrated navigation was verified for hypersonic velocity by the HWIL simulation of hypersonic vehicles [4], as shown in Figure 11. e HWIL simulation equipment includes the real-time simulator, GNC analysis computer, 3D display, three-axis rotation table, IMU simulator, GNSS simulator, product interface system, and actuator load simulator. e on-board equipment includes the on-board computer, SINS/BDS system, and actuator system.
We propose a classic boost-glide flight trajectory [22,23]. e equation describing this trajectory is detailed in [22], and the trajectory data are obtained from the simulation platform in [23]. e initial parameters of the   trajectory are shown in Table 1, and the simulation diagram of the trajectory is shown in Figure 12.

Conclusions
In this work, we considered the hypersonic boost-glide vehicles as the object, proposed the SINS mechanization and SINS numerical update algorithm in the LCI frame, and designed the SINS/BDS integrated algorithm in the LCI frame. In addition, we conducted vehicle field tests and hardware-in-the-loop simulations, and the accuracy of the simulation results was found to meet the navigation requirements. When the pitch angle is 90°, the attitude angle in the LCI frame is not singular, which meets the requirements for the vertical launch of the hypersonic boost-glide vehicles. erefore, the proposed integrated navigation system is especially suitable for hypersonic boost-glide vehicles and other vertically launched aircrafts.
Data Availability e data used in this article are not available for public disclosure.

Conflicts of Interest
e authors declare that they have no conflicts of interest. Mathematical Problems in Engineering 15