Information Fusion Based on Complementary Filter for SINS/CNS/GPS Integrated Navigation System of Aerospace Plane

In order to solve the problems of heavy computational load and poor real time of the information fusion method based on the federated Kalman filter (FKF), a novel information fusion method based on the complementary filter is proposed for strapdown inertial navigation (SINS)/celestial navigation system (CNS)/global positioning system (GPS) integrated navigation system of an aerospace plane. The complementary filters are designed to achieve the estimations of attitude, velocity, and position in the SINS/CNS/GPS integrated navigation system, respectively. The simulation results show that the proposed information fusion method can effectively realize SINS/CNS/GPS information fusion. Compared with FKF, the method based on complementary filter (CF) has the advantages of simplicity, small calculation, good real-time performance, good stability, no need for initial alignment, fast convergence, etc. Furthermore, the computational efficiency of CF is increased by 94.81%. Finally, the superiority of the proposed CF-based method is verified by both the semi-physical simulation and real-time system experiment.


Introduction
The aerospace plane is a reusable aerospace vehicle with the features of aircraft, spacecraft, and carrier, in which the accuracy and reliability requirements of the navigation system are very high due to its wide flight envelope, complex flight environment, and special mission [1]. The common navigation systems in aerospace vehicle include the strapdown inertial navigation system (SINS), global positioning system (GPS), celestial navigation system (CNS), and so on. Considering the black-out phenomenon in atmospheric reentry, the single SINS navigation using single-axis rotation modulation [2] or a second-order damping scheme [3] was adopted for the navigation of hypersonic vehicles. In order to serve the needs of spacecraft, the global navigation satellite system (GNSS) based on GPS [4], global navigation satellite system (GLONASS) [5], GPS/Galileo [6], or GPS/beidou navigation satellite system (BDS) [7] was proposed. An approach using the measurements from GNSS was proposed for the navigation of circumlunar spacecraft [8] or high-Earth-orbits spacecraft [9]. References [10,11] studied the high-precision autonomous celestial navigation based on astronomical information for spacecraft.
In view of the advantages and disadvantages of each navigation subsystem, it is difficult for any single navigation system to provide high-precision navigation parameters independently. In this background, the method of integrated navigation was proposed. Compared with a single navigation system, the integrated navigation system has higher accuracy as well as better reliability and information completeness; therefore, it can be applied for the scenario of an aerospace plane.

SINS Mechanization
The navigation parameters of SINS include the attitude quaternion q n b = q 0 q 1 q 2 q 3 T , velocity v n = v E v N v U T , and position p n = L λ h T ; the superscript n denotes the navigation frame and the subscript b denotes the body-fixed frame, q i (i = 0, 1, 2, 3) denotes the ith component of attitude quaternion. v E , v N , and v U represent the eastern velocity, northern velocity, and vertical velocity, respectively. L, λ, and h represent the latitude, longitude, and altitude, respectively, and the mechanization equations for SINS are expressed in Equation (1).
where ω b nb represents the angular velocity of the b frame relative to the n frame expressed in the b frame, f n represents the measurement of a specific force in the n frame, ω n ie represents the Earth rotation rate denoted in the n frame, ω n en represents the angular velocity of the n frame relative to the e frame coordinated in the n frame, g n represents the gravity acceleration expressed in the n frame, R M and R N are the radii of curvature in meridian and prime vertical, (L, λ, h) is the position of the vehicle in latitude, longitude and altitude.

CNS Measurement Transformation
The output of the CNS is the attitude quaternion q i,cns b from the b frame to the i frame, but the output of SINS is the attitude quaternion q n,ins b from the b frame to the n frame. If the complementary filter is selected for the information fusion of SINS/CNS integration, the parameter conversion is required firstly; that is to say, q i,cns b (the attitude quaternion from the b frame to the i frame) must be converted into q n,cns b (the attitude quaternion from the b frame to the n frame). The algorithm of CNS quaternion transformation from q i,cns b to q n,cns b is given in this subsection. According to the chain rule, the attitude quaternion q n,cns b is q n,cns b = q n e ⊗ q e i ⊗ q i,cns b (2) where q n e is the transformation quaternion from the e frame to the n frame, and q e i is the transformation quaternion from the i frame to the e frame.
cos(λ gps ) 0 − sin(L gps ) cos(λ gps ) − sin(L gps ) sin(λ gps ) cos(L gps ) cos(L gps ) cos(λ gps ) cos(L gps ) sin(λ gps ) sin(L gps ) According to the relationship between the direction cosine matrix and corresponding quaternion, the expression of quaternion q n e is as follows: q n e = Mat2Quat(C n,gps e ) (4) where Mat2Quat(•) represents the function that calculates the quaternion according to the direction cosine matrix, and its expression is shown in Equation (5).
For any direction cosine matrix C = T ij 3×3 (i, j = 1, 2, 3), the corresponding attitude quaternion is defined as q q 0 q 1 q 2 q 3 T ; then, the method of calculating q i (i = 0, 1, 2, 3) is as follows: where sign(•) represents the sign function.

The Method for Calculating q e i
The quaternion q e i can be calculated according to the Greenwich Mean Time and Earth rotation rate. Given Greenwich Mean Time t G and Earth rotation rate ω ie , the Greenwich Hour Angle G ST (the rotation angle of the e frame relative to the i frame along the z i axis) is obtained by multiplying ω ie and t G together: Then, the transformation matrix C e i can be calculated, as shown in Equation (7): Sensors 2020, 20, 7193 5 of 26 According to the relationship between the direction cosine matrix and corresponding attitude quaternion, the quaternion q e i can be calculated using C e i , and its expression is shown in Equation (8):

The Principle of Complementary Filter
The complementary filtering method can restrain the error divergence of the integrated navigation system, which is designed according to the frequency-domain characteristics of errors of different subsystems. For any signal variable X, when it is measured using two uncorrelated methods (for convenience, they are denoted as Method #1 and Method #2, respectively), two measurement results Z 1 and Z 2 can be obtained, and they can be expressed as: where Z 1 denotes the measurement result using Method #1, V 1 denotes the measuring error of Z 1 ; Z 2 denotes the measurement result using Method #2, V 2 denotes the measuring error of Z 2 ; V 1 and V 2 are assumed to be the low-frequency error and high-frequency error, respectively. In order to obtain an accurate estimate of variable X, the measurements Z 1 and Z 2 are filtered using the high-pass filter (H.P.F. for short) and low-pass filter (L.P.F for short), respectively, so that the errors V 1 and V 2 could be eliminated.
In the design of a complementary filter, the transfer functions of H.P.F. and L.P.F are denoted as G 1 (s) and G 2 (s), respectively. By selecting appropriate parameters for H.P.F. and L.P.F, the condition G 1 (s) + G 2 (s) = 1 can be satisfied; then, the estimateX(s) of X(s) is obtained by applying the signal reconstruction (S.R.) method (i.e., the sum of the outputs of G 1 (s) and G 2 (s)). For the low-pass filter G 2 (s) and high-pass filter G 1 (s), the appropriate cut-off frequencies are selected so that G 1 (s)V 1 (s) and G 2 (s)V 2 (s) are approximate to zero, and the sum of the outputs of G 1 (s) and G 2 (s) will be close to the true value of signal variable X, as shown in Equation (10): The principle of the complementary filter is shown in Figure 1.

Design of Complementary Filters for SINS/CNS/GPS Integration
There are two sources of attitude information for SINS/CNS/GPS integrated navigation: the attitude in the n frame calculated by SINS (denoted as q n,ins b ), the attitude in the i frame measured by CNS (denoted as q i,cns b ), and the velocity information sources include the velocity in the n frame calculated by SINS (denoted as v n ins ) and the velocity in the n frame measured by GPS (denoted as v n gps (s)). Similarly, the position information sources include the position in the n frame calculated by SINS (denoted as p ins ) and the position in the n frame measured by GPS (denoted as p gps ). The navigation information (including the attitude, velocity, and position information) obtained by SINS has the advantages of good dynamic performance, high real-time stability, and good short-term stability, but the navigation errors accumulate over time, and the long-term stability is poor. The error sources of SINS mainly include the gyro drifts/accelerometer biases, iterative algorithm errors, initial alignment errors, etc. Due to the cumulative effect of integral calculation, the errors of SINS navigation parameters, which include the attitude error, velocity error, and position error, will accumulate over time, and all the frequency spectrums are mainly in the low frequency segment. The navigation information measured by CNS and GPS (including the attitude from CNS, the velocity and position from GPS) has the advantages of high precision, error convergence, and good long-term stability, but the data update rate is low, and the real-time performance is poor. The measurement errors of CNS and GPS do not accumulate over time and can be approximated as Gaussian white noise. By analyzing the frequency-domain characteristics of the navigation errors of each subsystem, the complementary filters are designed to realize the information fusion of an SINS/CNS/GPS integrated navigation system, and the navigation results with higher precision are obtained. The process of complementary filtering information fusion for an SINS/CNS/GPS integrated navigation system is shown in Figure 2.

Design of Complementary Filters for SINS/CNS/GPS Integration
There are two sources of attitude information for SINS/CNS/GPS integrated navigation: the attitude in the n frame calculated by SINS (denoted as (denoted as gps p ). The navigation information (including the attitude, velocity, and position information) obtained by SINS has the advantages of good dynamic performance, high real-time stability, and good short-term stability, but the navigation errors accumulate over time, and the long-term stability is poor. The error sources of SINS mainly include the gyro drifts/accelerometer biases, iterative algorithm errors, initial alignment errors, etc. Due to the cumulative effect of integral calculation, the errors of SINS navigation parameters, which include the attitude error, velocity error, and position error, will accumulate over time, and all the frequency spectrums are mainly in the low frequency segment. The navigation information measured by CNS and GPS (including the attitude from CNS, the velocity and position from GPS) has the advantages of high precision, error convergence, and good long-term stability, but the data update rate is low, and the real-time performance is poor. The measurement errors of CNS and GPS do not accumulate over time and can be approximated as Gaussian white noise. By analyzing the frequency-domain characteristics of the navigation errors of each subsystem, the complementary filters are designed to realize the information fusion of an SINS/CNS/GPS integrated navigation system, and the navigation results with higher precision are obtained. The process of complementary filtering information fusion for an SINS/CNS/GPS integrated navigation system is shown in Figure 2.

Complementary Filter for Attitude Estimation
As shown in Section 2.3, the attitude calculated by SINS is  , and the attitude measured by CNS is q i,cns b . By using the position information of GPS, the attitude of CNS in the n frame can be calculated, as shown in Equation (11): where q n,gps e represents the transformation quaternion from the e frame to the n frame; q e i represents the transformation quaternion from the i frame to the e frame; q n,gps e can be calculated based on the position of GPS; q e i can be calculated based on the UTC time and EOP parameters, in which the EOP parameters represent the Earth orientation parameters provided by the International Earth Rotation and Reference Systems Service (IERS).
For the attitude estimation, the complementary filter is denoted as CF1, which is composed of a first-order high-pass filter and a first-order low-pass filter. The transfer function of a first-order high-pass filter is denoted as G H 1 (s) = s/ s + ω c 1 , and the transfer function of a first-order low-pass filter is denoted as G L 1 (s) = ω c 1 / s + ω c 1 , where ω c 1 is the cut-off angular frequency of CF1, and the cut-off frequency corresponding to ω c 1 is f c 1 = ω c 1 /(2π). Then, the attitude estimation in the s-domain is as follows: Suppose that the filtering period is T f = t k − t k−1 . Discretize Equation (12) to obtain the recursive form in time domain for attitude estimation, as shown in Equation (13): where ω c 1 represents the angular frequency corresponding to CF1 s cut-off frequency f c 1 = ω c 1 /(2π).

Complementary Filter for Velocity Estimation
As mentioned, the velocity of SINS is v n ins , and the velocity of GPS is v n gps . For the velocity estimation, the complementary filter is denoted as CF2, which is composed of a first-order high-pass filter and a first-order low-pass filter. The transfer function of a first-order high-pass filter is denoted as G H 2 (s) = s/ s + ω c 2 , and the transfer function of a first-order low-pass filter is denoted as G L 2 (s) = ω c 2 / s + ω c 2 , where ω c 2 represents the angular frequency corresponding to CF2 s cut-off frequency f c 2 = ω c 2 /(2π). Then, the velocity estimation in the s-domain is as follows: Similarly, the filtering period is also T f . Discretize Equation (14) to obtain the recursive form in the time domain for velocity estimation, as shown in Equation (15):

Complementary Filter for Position Estimation
For position estimation, the position of SINS is p ins , the position of GPS is p gps , and the complementary filter is denoted as CF3, which is also composed of a first-order high-pass filter and a first-order low-pass filter. The transfer function of a first-order high-pass filter is G H 3 (s) = s/ s + ω c 3 , and the transfer function of a first-order low-pass filter is G L 3 (s) = ω c 3 / s + ω c 3 , where ω c 3 is the angular frequency corresponding to CF3 s cut-off frequency f c 3 = ω c 3 /(2π). Then, the position estimation in the s-domain is as follows:p (s) = G H 3 (s)p ins (s) + G L 3 (s)p gps (s). (16) Similarly, take the filtering period T f as a discrete step, and discretize Equation (16) to obtain the recursive form in the time domain for position estimation, as shown in Equation (17): Each time after the process of complementary filtering is completed, the navigation parameters of SINS are corrected by feedback of the estimation results, as shown in Equation (18):

Error Analysis in Frequency Domain and Selection of Cut-Off Frequency
The cut-off frequency is the most key parameter of a complementary filter, which determines the effect of information fusion. In order to obtain good estimation results, it is necessary to select the cut-off frequencies properly for the complementary filters by analyzing the errors of each navigation subsystem in the frequency domain. According to the principle of complementary filter in Section 3.1, the analysis shows that the selection of cut-off frequency of the complementary filter determines the influential weight of each navigation subsystem (SINS, CNS, or GPS) on the SINS/CNS/GPS integrated navigation system. In conclusion, the larger the cut-off frequency f c i (i = 1, 2, 3), the greater the influence of the measurement noise of the auxiliary subsystem (CNS or GPS) on the estimation result, the smaller the influence of the errors of SINS, and vice versa. During the on-orbit flight of an aerospace plane, the navigation errors of CNS and GPS are mainly the measurement noises, which can be approximated as Gaussian white noise. Under the premise of satisfying the dynamic performance requirements of SINS, the cut-off frequency of the complementary filter should be as small as possible, so that the measurement noise errors of CNS and GPS can be filtered out in the largest frequency range. In this paper, the errors of each subsystem are analyzed in the time and frequency domain. According to the relation between the cut-off frequencies and statistical errors, the cut-off frequencies of complementary filters for attitude/velocity/position estimation are optimally selected, respectively.

Selection of Cut-Off Frequency Based on Nonlinear Optimization Theory
In this subsection, the analysis method based on power spectral density is used to analyze the influence of the noise of each subsystem on the output noise of the SINS/CNS/GPS integrated navigation system theoretically, and an optimal selection scheme of the cut-off frequency of the complementary filter based on the constrained nonlinear optimization method is proposed, so as to provide a theoretical reference for the selection of the optimal cut-off frequency.
In SINS/CNS/GPS integrated navigation, the measurement noises of navigation sensors such as gyro, accelerometer, CNS, and GPS are all stationary random noise. For the convenience of theoretical analysis, these noises are assumed to be Gaussian white noise with zero mean. In this paper, the measurement noises of a gyro, accelerometer, CNS, GPS velocity, and GPS position are denoted as w gro (t), w acc (t), w cns_att (t), w gps_vel (t), and w gps_pos (t), respectively. The random walk noises in SINS attitude, velocity, and position are denoted as w ins_att (t), w ins_vel (t), and w ins_pos (t). q w gro , q w acc , q w cns_att , q w gps_vel , and q w gps_pos represent the variance intensities of gyro angular velocity measurement noise, accelerometer specific force measurement noise, CNS attitude measurement noise, GPS velocity measurement noise, and GPS position measurement noise, respectively. δ() represents the Dirac delta Sensors 2020, 20, 7193 9 of 26 function. Then, the covariances and autocorrelation functions of measurement noises of navigation sensors are expressed as: where µ = t − τ represents the time interval between the time instants t and τ.
The power spectral density functions of measurement noises of these navigation sensors are where F{} is the Fourier transform.
For the convenience of analysis, the motion state of the vehicle and the random error propagation model of SINS are simplified. Assuming that the vehicle is in an on-orbit cruise phase and its motion state is uniform straight-line flight, then the error model of SINS can be simplified as a linear time-invariant system.
The attitude random walk noise w ins_att (t) of SINS is the integral of gyro measurement noise w gro (t), and the frequency response function of the integrator is 1/( jω). According to the response characteristics of the linear system to stationary process, the power spectral density of w ins_att is Similarly, the SINS velocity random walk noise w ins_vel is the integral of measurement noise of the accelerometer, and the power spectral density of w ins_vel is Due to the measurement noises of inertial sensors and the calculation error of the SINS update algorithm, there will be white noise in the velocity information of SINS, which could be denoted as w ins_nv , and its variance intensity is denoted as q w ins_nv . The power spectral density of w ins_nv is S w ins_nv (ω) = q w ins_nv , and the position random walk noise w ins_pos of SINS is the integral of velocity white noise; its power spectral density is as follows: In this paper, the random errors of CF1, CF2, and CF3 are denoted as w CF_att , w CF_vel , and w CF_pos , respectively. Firstly, the power spectral density of a random error of attitude estimation complementary filter CF1 is analyzed. According to the design of complementary filters in Section 3.2, the transfer functions of the high-pass filter and low-pass filter of CF1 are rewritten as follows: From the transfer functions, the corresponding frequency response functions can be obtained as follows: According to the principle of complementary filter, the random error of CF1 is the sum of two random errors, of which one is the output of G H 1 ( jω) under the input of SINS attitude random walk noise, and the other is the output of G L 1 ( jω) under the input of CNS attitude measurement noise. Since SINS and CNS are two independent subsystems, the random error of SINS is orthogonal to that of CNS, so the power spectral density S w CF_att (ω) of the random error of CF1 is (26), we can obtain By substituting ω = 2π f into Equation (27), the expression of the power spectral density function of w CF_att with respect to frequency f can be obtained.
In a SINS/CNS/GPS integrated navigation system, the maximum frequency of a filter's bandwidth is limited due to the discrete sampling of a signal. Among the three subsystems of SINS, CNS, and GPS, the sampling frequency f s ins of SINS is the highest. According to the sampling theorem of a stationary random process, the power spectral density of the output noise of a SINS/CNS/GPS integrated navigation system only needs to be considered in the frequency range of 0, f s ins /2 . By calculating the definite integral of S w CF_att ( f ) with respect to f from 0 to f ins s /2, the average power P w CF_att of the random error w CF_att of CF1 in the frequency band of 0, f ins s /2 could be obtained, as shown in Equation (29): By substituting ω c 1 = 2π f c 1 into Equation (29), the expression of average power P w CF_att with respect to the cut-off frequency f c 1 can be obtained, as shown in Equation (30): In the algorithm design of a practical integrated navigation system, the parameters of f ins s , q w gro , and q w cns_att are constant, and P w CF_att is a nonlinear function of f c 1 . The optimal value of f c 1 is equivalent to the value that minimizes P w CF_att . By differentiating P w CF_att with respect to f c 1 , we can obtain It can be seen from Equation (31) that the derivative of P w CF_att with respect to f c 1 has a very complex form, so it is difficult to find the minimum value of P w CF_att in analytical form by using the stationary point method. Considering that in a specific SINS/CNS/GPS integrated navigation system, the values of f ins s , q w gro , q w cns_att are known or can be obtained by calibration, thus, the minimum value of P w CF_att with respect to f c 1 can be solved by the optimization method, and the nonlinear optimization function of the MATLAB Optimization Toolbox can be used to obtain the numerical solution of the optimization problem.
Combined with optimization theory, the optimal solution of f c 1 can be regarded as a constrained nonlinear optimization problem, and the optimal mathematical model is established as where f c L is the infimum of the range of the CF's cut-off frequency, which can be determined by analyzing the error propagation characteristics of the SINS. Under the condition of f c 1 ≥ f c L , the high-pass filter of CF1 can effectively restrain the oscillation and divergence low-frequency errors in the SINS attitude. Many studies about inertial navigation analyzed the error propagation characteristics of the inertial navigation system, which are not outlined here. Furthermore, the selection of f c 1 should meet the requirement of "minimizing the influence of measurement noise of CNS", thus f c 1 ≤ f cns s /2, i.e., the constraint condition in Equation (32). For programming implementation, the MATLAB built-in function "fmincon" can be used to solve Equation (32).
The power spectral density analyses of random errors of CF2 and CF3 are similar to that of CF1. The detailed analysis process will not be introduced here, the analysis results of random noises of CF2 and CF3 are given directly (including the power spectral density, average power, and optimal mathematical model).
According to the power spectral density analysis of random noise of CF2, the following conclusions can be obtained: Power spectral density of random error of CF2 S w CF_vel ( f ): Average power of random error of CF2 P w CF_vel : In order to solve the optimization problem of f c 2 , the mathematical model of constrained nonlinear optimization is established as follows: where f c L is the spectrum bandwidth of INS errors. Under the condition of f c 2 ≥ f c L , the high-pass filter of CF2 can effectively eliminate the oscillation and divergence terms of SINS velocity errors. f gps s is the sampling frequency of GPS.
According to the power spectral density analysis of random noise of CF3, the following conclusions can be obtained: Power spectral density of random error of CF3 S w CF_pos ( f ): Sensors 2020, 20, 7193

of 26
Average power of random error of CF3 P w CF_pos : To solve the optimization problem of f c 3 , the mathematical model of constrained nonlinear optimization is established as follows: where f c L is the spectrum bandwidth of SINS errors. Under the condition of f c 3 ≥ f c L , the high-pass filter of CF3 can effectively eliminate the oscillation and divergence terms in INS position errors. f gps s is the sampling frequency of GPS.

Spectrum Analysis Based on DFT
According to the theory of digital signal processing, the discrete Fourier transform (DFT) is often used for signal spectrum analysis. If an N-point time-domain sampling sequence is denoted as x(n) (0 ≤ n ≤ N − 1), then the corresponding N-point DFT of x(n) is defined as: In order to study the frequency distribution of error, DFT is often used to analyze the frequency spectrum of an error signal, mainly focusing on the amplitude-frequency characteristics. For DFT, the sampling frequency of the signal is denoted as f s ; then, the steps of Fourier analysis on the N-point sampling data x(n) are as follows: Step 1: The DFT transform on the N-point sampling data x(n) is performed to obtain the N-point DFT X(k).
Step 2: Calculate the spectral amplitude X(k) for the spectrum data X(k) of the ordinal number 0 ∼ f ix(N/2), where f ix(•) denotes the function rounding a number to the nearest integer.
Step 3: Divide all amplitude data of the spectra obtained in Step 2 by N; then, multiply the amplitude data of ordinal number 2 ∼ f ix(N/2) + 1 by 2, convert the two-sided spectrum into a one-sided spectrum, and plot the amplitude-frequency curve.
Step 4: The frequencies of spectral lines labeled on the abscissa axis of the amplitude-frequency figure are sequentially k f s /N Hz, ( k = 0 ∼ f ix(N/2)).
Step 5: The ordinate of spectral line corresponds to the harmonic amplitude of sampling signal x(n) at the corresponding frequency, the unit of amplitude is the same as that of time-domain sampling data, and the amplitude of DC component is at zero frequency.

Selection of Cut-Off Frequencies Based on Simulation Test
In order to determine the cut-off frequencies of complementary filters for the attitude/velocity/position estimation in an SINS/CNS/GPS integrated navigation system, the frequency-domain analyses were conducted on the measurement errors of SINS, CNS, and GPS subsystems. The analysis results, which include the error curves, amplitude-frequency diagrams, and relationship graphs between the estimation errors and selected cut-off frequencies, are shown in Figures 3-11, respectively.
Firstly, analyze and select the cut-off frequency of the complementary filter for attitude estimation (i.e., the cut-off frequency of CF1). The attitude errors of SINS and CNS subsystems obtained in the static test are analyzed in time and frequency domain, and the analysis results are shown in Figures 3-5 where φ x , φ y , and φ z denote the pitch, roll, and yaw misalignment angles, respectively, the unit of which is arcsecond ("); t denotes the time axis, and its unit is second (s); φ x , φ y and φ z denote the amplitudes of the pitch, roll, and yaw error spectrums, the unit of which is arcsecond ("); f denotes the frequency axis, its unit is Hertz (Hz); σ(φ x ), σ φ y , and σ(φ z ) denote the standard deviations of attitude errors in SINS/CNS integration, and f c 1 denotes the cut-off frequency of CF1. Figure 3 shows the time curves and amplitude spectrums of CNS attitude errors, of which the three subfigures on the left are the time curves, and three subfigures on the right are the amplitude spectrums. Figure 4 shows the time curves and amplitude spectrums of SINS attitude errors. Similarly, the three subfigures on the left are the time curves, and the three subfigures on the right are the amplitude spectrums. It can be seen from Figure 3 to Figure 4 that the attitude errors of CNS do not diverge and are less than 100", while the time curves of attitude errors are characterized by disorderly jumping noise. According to the corresponding amplitude spectrums, the errors of CNS can be approximated as white noise. The attitude errors of SINS have a tendency to oscillate and diverge, which are close to 20 in 3600 s. From the amplitude spectrums, it can be seen that the attitude errors of SINS are mainly distributed in the low frequency band of 0-1 Hz. In order to obtain the statistical relationship between the cut-off frequency f c 1 and attitude estimation accuracy σ(φ i ), let f c 1 take different values, and Monte-Carlo simulation is performed to record the corresponding values of σ(φ i ) (i = x, y, z) and Hz. According to the statistical results of Monte-Carlo simulation, the relationships between f c 1 and σ(φ i ) (i = x, y, z) are obtained. As shown in Figure 5, the effect of cut-off frequency on the attitude error can be seen. The data pairs including the minimum values of σ(φ i ) (i = x, y, z) and corresponding f c 1 are pitch (0.09 Hz, 0.2116'), roll (0.09 Hz, 0.2140'), and yaw (0.12 Hz, 0.2121'), which are marked with red asterisks(*) in Figure 5. It can be seen that the cut-off frequencies corresponding to the minimum values of three attitude errors are very close. In order to reduce the algorithmic complexity and computational load, the same cut-off frequency values are selected for the pitch, roll, and yaw estimations, i.e., f c Sensors 2020, 20, x FOR PEER REVIEW 14 of 32 x φ , y φ and z φ denote the amplitudes of the pitch, roll, and yaw error spectrums, the unit of which is arcsecond (″); f denotes the frequency axis, its unit is Hertz ( Hz ); ( )      In order to select the appropriate cut-off frequency of a complementary filter for velocity estimation (i.e., the cut-off frequency of CF2, f c 2 ), the velocity errors of the SINS and GPS subsystems obtained in the static test are analyzed in the time and frequency domains, and the analysis results are shown in Figures 6-8, where δV E , δV N , and δV U denote the velocity errors in the East, North, and Up directions, respectively, the unit of which is meters per second (m/s); t denotes the time axis, and its unit is second (s); |δV E |, |δV N |, and |δV U | denote the amplitude spectrums of velocity errors in the East, North, and Up directions, the unit of which is meters per second (m/s); f denotes the frequency axis, its unit is Hertz (Hz); σ(δV E ), σ(δV N ), and σ(δV U ) denote the standard deviations of velocity errors in SINS/GPS integration, and f c 2 denotes the cut-off frequency of CF2. Figure 6 shows the time curves and amplitude spectrums of GPS velocity errors, of which three subfigures on the left are the time curves, and three subfigures on the right are the amplitude spectrums. Figure 7 shows the time curves and amplitude spectrums of SINS velocity errors. Similarly, three subfigures on the left are the time curves, and three subfigures on the right are the amplitude spectrums. It can be seen from Figures 6 and 7 that the velocity errors of GPS do not diverge and are all less than 0.5 m/s, the time curves of GPS velocity errors are characterized by disorderly jumping noise. According to the corresponding amplitude spectrums, the velocity errors of GPS can be approximated as Gaussian white noise. The velocity errors of SINS have a tendency to oscillate and diverge, which are close to 20 m/s in 3600 s. From the amplitude spectrum, it can be seen that the velocity errors of SINS are mainly distributed in the low frequency band of 0 to 0.5 Hz. In order to obtain the statistical relationship between the cut-off frequency value and velocity estimation accuracy, let the cut-off frequency f c 2 take different values, and Monte-Carlo simulation is performed to record the values of σ(δV i ) (i = E, N, U) and corresponding f c 2 . For Monte-Carlo simulation of CF2, the selection of f c 2 value points is the same as that of f c 1 . According to the statistical results of Monte-Carlo simulation, the relationships between f c 2 and σ(δV i ) (i = E, N, U) are obtained, as shown in Figure 8, so the effect of cut-off frequency on the error of velocity estimation can be seen. The data pairs including the minimum values of σ(δV i ) and corresponding f c 2 are East (0.03 Hz, 0.0413m/s), North (0.03 Hz, 0.0408m/s), and Up (0.03 Hz, 0.0381m/s), which are marked with red asterisks(*) in Figure 8. It can be seen that the cut-off frequencies corresponding to the minimum values of the three velocity errors are very close. Therefore, the same cut-off frequency value can be selected for velocity estimation; there, it is selected as f c 2 = 0.03 Hz, ω c 2 = 2π f c 2 = 0.06π.   Figure 9 shows the time curves and amplitude spectrums of GPS position errors, of which three subfigures on the left are the time curves, and three subfigures on the right are the amplitude spectrums. Figure  Finally, the cut-off frequency is selected for position estimation (i.e., the cut-off frequency of CF3, f c 3 ). The position errors of SINS and GPS subsystems obtained in the static test are analyzed in time and frequency domain, and the analysis results are shown in Figures 9-11. Thereinto, δL, δλ, and δH denote the latitude, longitude, and altitude errors, respectively, the unit of which is meter (m); t denotes the time axis, and its unit is second (s); |δL|, |δλ|, and |δH| denote the amplitude spectrums of latitude, longitude, and altitude errors, the unit of which is meter (m); f denotes the frequency axis, its unit is Hertz (Hz); σ(δL), σ(δλ), and σ(δH) denote the standard deviations of position errors in SINS/GPS integration; and f c 3 denotes the cut-off frequency of CF3. Figure 9 shows the time curves and amplitude spectrums of GPS position errors, of which three subfigures on the left are the time curves, and three subfigures on the right are the amplitude spectrums. Figure 10 shows the time curves and amplitude spectrums of SINS position errors. Similarly, the three subfigures on the left are the time curves, and the three subfigures on the right are the amplitude spectrums. It can be seen from Figures 9 and 10 that the position errors of GPS do not diverge and are less than 30 m, while the time curves of the GPS position errors are characterized by disorderly jumping noise. According to the corresponding amplitude spectrums, the position errors of GPS can be approximated as Gaussian white noise. The position errors of SINS will oscillate and diverge, which are close to the order of magnitude of 2 × 10 4 m at 3600 s. As seen from the amplitude spectrum, the position errors of SINS are mainly distributed in the low frequency band of 0 to 2 Hz. Similarly, in order to obtain a statistical relationship between the cut-off frequency and position estimation accuracy, let the cut-off frequency f c 3 take different values, and Monte-Carlo simulation is performed to record the values of σ(δP i ) (i = L, λ, H) and corresponding f c 3 ; the selection of f c 3 value points is the same as that of f c 1 . According to the statistical results of Monte-Carlo simulation, the relationships between f c 3 and σ(δP i ) (i = L, λ, H) are obtained. As shown in Figure 11, the effect of cut-off frequency on the position errors of SINS/GPS integration can be seen. The data pairs including the minimum values of σ(δP i ) and corresponding f c 3 are latitude (0.15 Hz, 6.7436m), longitude (0.15 Hz, 5.6282m), and altitude (0.09 Hz, 6.2559m), which are marked with red asterisks(*) in Figure 11. The cut-off frequencies corresponding to the minimum values of the latitude, longitude, and altitude errors are very close. Based on the analysis above, the value of f c 3 can be selected as f c 3 = 0.12 Hz, ω c 3 = 2π f c 3 = 0.24π.
Sensors 2020, 20, x FOR PEER REVIEW 20 of 32 marked with red asterisks(*) in Figure 11. The cut-off frequencies corresponding to the minimum values of the latitude, longitude, and altitude errors are very close. Based on the analysis above, the The minimum standard deviations of navigation parameters and corresponding cut-off frequencies are summarized in Table 1. Finally, the values of f c 1 , f c 2 , and f c 3 are determined as 0.1, 0.03, and 0.12 Hz, respectively. Table 1. Statistics of minimum standard deviations and corresponding cut-off frequencies.
Errors  Figure 11. The relationships between f c 3 and σ(δP i ).

Simulation and Experiment
The real-time data of an SINS/CNS/GPS integrated navigation system of an aerospace plane in flight are difficult to obtain. In this paper, the method of trajectory simulation combined with the error characteristics of a real sensor's outputs is used to generate the simulation data of each navigation subsystem, and the semi-physical simulation of an SINS/CNS/GPS integrated navigation algorithm based on CF is performed to verify the validity. The semi-physical simulation experiment using this system has important practical significance for studying the characteristics of a SINS/CNS/GPS integrated navigation system of spacecraft under actual noise. In order to fully verify the priority of CF in real-time performance, an experimental test based on a real-time system was done.

Design and Implementation of Semi-Physical Simulation System
Combined with the application background of spacecraft, the semi-physical simulation platform of the SINS/CNS/GPS integrated navigation system is constructed, which consists of the hardware part (including IMU, GPS, CNS, navigation computer, etc.) and software part (including the trajectory simulation, integrated navigation algorithm, SINS updating algorithm, etc.). By modeling the sampling data of SINS, GPS, and CNS, the error characteristics of each navigation subsystem are obtained. According to the presupposed maneuver conditions, the required error-free data of navigation sensors are generated by using the method of trajectory simulation. According to the error characteristics of each subsystem, the real outputs of subsystems can be simulated by adding errors to the corresponding error-free data. Finally, the algorithm software of SINS and SINS/CNS/GPS integration are performed, and the navigation results are output.
The flow chart of semi-physical simulation for an SINS/CNS/GPS integrated navigation system is shown in Figure 12. The process is as follows: (1) Set the initial parameters of trajectory simulation, perform the trajectory simulation, and generate the error-free trajectory data (including the error-free attitude, velocity, and position of spacecraft). (2) According to the obtained trajectory data, calculate the error-free output of IMU. Collect the IMU data in a static test, and subtract the corresponding mean value from the IMU data to obtain the error data of inertial sensors. Then, add the error data to the error-free IMU data, thus obtaining the simulated output data of IMU with real errors. Then, calculate the attitude, velocity, and position of the SINS. (3) The speed and position data of GPS in the static test are collected, and the means are subtracted from collected data to obtain the error data of GPS. Add the GPS's error data to the error-free position/velocity data of spacecraft from trajectory simulation generated in step (1), thereby obtaining the simulated output data of GPS with real errors. (4) Similarly, the error data of the CNS can be obtained by subtracting the true value of attitude from the collected off-line attitude data of CNS. Then, add the error data of the CNS to the error-free attitude data of the spacecraft from the trajectory simulation generated in step (1), thereby obtaining the simulated output data of CNS with real errors. (5) By using the data from navigation subsystems generated in steps (2) (4) Similarly, the error data of the CNS can be obtained by subtracting the true value of attitude from the collected off-line attitude data of CNS. Then, add the error data of the CNS to the error-free attitude data of the spacecraft from the trajectory simulation generated in step (1), thereby obtaining the simulated output data of CNS with real errors.
(5) By using the data from navigation subsystems generated in steps (2), (3), and (4), the information fusion based on complementary filter is performed to obtain the optimal estimations of attitude/velocity/ position parameters of an SINS/CNS/GPS integrated navigation system.

Simulation Results
The semi-physical simulation of the proposed CF-based SINS/CNS/GPS integrated navigation algorithm is performed, and the navigation errors obtained by comparing the integrated navigation results with the corresponding trajectory parameters generated by the trajectory generator are shown in  In order to further analyze the performance of the proposed algorithm, the semi-physical simulation of SINS/CNS/GPS integrated navigation algorithm based on FKF (federated Kalman filter) is carried out under the same simulation conditions. The estimation errors of CF are compared with that of FKF, as shown in Table 2. For both CF and FKF, the mean value of each navigation parameter error is one order of magnitude smaller than the standard deviation of the corresponding error. The means of errors are no longer listed in the table, and only the standard deviations of errors are compared and analyzed. The simulation results show that the statistical errors of the CF-based integrated navigation algorithm are close to that of the FKF-based integrated navigation algorithm. In terms of filtering accuracy, the estimation errors of the two filtering methods are in the same order of magnitude; however, in terms of the real-time and complexity, because FKF needs to perform a large number of high-order matrix operations, the computational complexity of FKF is much higher than that of CF, and the real-time performance of FKF is much poorer than that of CF. From Figure 13, Figure 14, Figure 15, when the CF processes of INS/CNS/GPS integration navigation start, the curves of attitude errors, velocity errors, and position errors converge quickly, and the noise components of navigation errors are basically filtered out. Throughout the process of simulation, the navigation parameters output by the CF-based SINS/CNS/GPS integrated navigation system remain convergent. The attitude errors are less than 0.7 , the velocity errors are less than 0.15 m/s, and the position errors are less than 10 m. The errors of navigation parameters are all within the acceptable range. So, the proposed method of CF-based information fusion can meet the requirements of long-term and high-precision navigation for an SINS/CNS/GPS integrated navigation system of an aerospace plane.
In order to further analyze the performance of the proposed algorithm, the semi-physical simulation of SINS/CNS/GPS integrated navigation algorithm based on FKF (federated Kalman filter) is carried out under the same simulation conditions. The estimation errors of CF are compared with that of FKF, as shown in Table 2. For both CF and FKF, the mean value of each navigation parameter error is one order of magnitude smaller than the standard deviation of the corresponding error. The means of errors are no longer listed in the table, and only the standard deviations of errors are compared and analyzed. The simulation results show that the statistical errors of the CF-based integrated navigation algorithm are close to that of the FKF-based integrated navigation algorithm. In terms of filtering accuracy, the estimation errors of the two filtering methods are in the same order of magnitude; however, in terms of the real-time and complexity, because FKF needs to perform a large number of high-order matrix operations, the computational complexity of FKF is much higher than that of CF, and the real-time performance of FKF is much poorer than that of CF.

Analysis of Real-Time Performance
In order to investigate the real-time performance of the proposed CF-based information fusion algorithm for SINS/CNS/GPS integrated navigation, the average single-time elapsed times of the FKF and CF are recorded by using TIC and TOC functions of MATLAB. Denote the average single-time elapsed times of the FKF and CF as t FKF and t CF , respectively. According to the recorded time, we have t CF = 3.3836 × 10 −5 s and t FKF = 6.5180 × 10 −4 s. In order to eliminate the effect resulted from the difference of computer performances, the relative elapsed times of FKF and CF are defined as: Obviously, t R FKF = 1, t R CF = 0.0519. The relative elapsed times of FKF and CF are shown in Figure 16. From Figure 16, it can be seen that for the information fusion of SINS/CNS/GPS integrated navigation, the elapsed time of the CF is much shorter than that of the FKF, and the average single-time elapsed time of the CF is just 5.19% of the average single-time elapsed time of the FKF. The above simulation results demonstrate that the proposed CF-based method of information fusion for multi-sensor integrated navigation system successfully overcomes the drawback of heavy computational load of the FKF-based method. When used for SINS/CNS/GPS integrated navigation, the proposed CF-based method can save 94.81% computational time compared with the FKF-based method.

Experimental Test Based on Real-Time System
In order to further verify the priorities of a proposed CF-based information fusion algorithm in terms of computational complexity and real-time performance, an experimental test based on a real-time system was conducted. The data processing architecture of a real-time system is shown in Figure 17, in which the navigation computer mainly includes Digital Signal Processor (DSP, Texas Instruments, TMS320C6748), Field Programmable Gate Array (FPGA, Altera, Cyclone V 5CEFA9F23I7N), Analog-to-Digital Converter (ADC, Analog Devices Inc., Norwood, MA, USA, AD9288), etc. The DSP chip is used to calculate the information fusion algorithms of SINS/CNS/GPS integrated navigation. The experimental process is as follows: Firstly, a frame of experimental data (including the data of SINS, CNS, and GPS) generated by semi-physical simulation is input into the DSP chip through the RS232/422 serial communication interface of the navigation computer. Subsequently, the two information fusion algorithms of SINS/CNS/GPS integrated navigation based on CF and FKF are carried out by the DSP chip, respectively. Eventually, the recorded time consumptions are output to a terminal display device.
recorded time, we have . In order to eliminate the effect resulted from the difference of computer performances, the relative elapsed times of FKF and CF are defined as: Obviously, The relative elapsed times of FKF and CF are shown in Figure 16. From Figure 16, it can be seen that for the information fusion of SINS/CNS/GPS integrated navigation, the elapsed time of the CF is much shorter than that of the FKF, and the average single-time elapsed time of the CF is just 5.19% of the average single-time elapsed time of the FKF.
The above simulation results demonstrate that the proposed CF-based method of information fusion for multi-sensor integrated navigation system successfully overcomes the drawback of heavy computational load of the FKF-based method. When used for SINS/CNS/GPS integrated navigation, the proposed CF-based method can save 94.81% computational time compared with the FKF-based method.

Experimental Test Based on Real-Time System
In order to further verify the priorities of a proposed CF-based information fusion algorithm in terms of computational complexity and real-time performance, an experimental test based on a real-time system was conducted. The data processing architecture of a real-time system is shown in Figure 17, in which the navigation computer mainly includes Digital Signal Processor (DSP, Texas Instruments, TMS320C6748), Field Programmable Gate Array (FPGA, Altera, Cyclone V 5CEFA9F23I7N), Analog-to-Digital Converter (ADC , Analog Devices Inc., Norwood, MA, USA, AD9288), etc. The DSP chip is used to calculate the information fusion algorithms of SINS/CNS/GPS integrated navigation. The experimental process is as follows: Firstly, a frame of experimental data (including the data of SINS, CNS, and GPS) generated by semi-physical simulation is input into the DSP chip through the RS232/422 serial communication interface of the navigation computer. Subsequently, the two information fusion algorithms of SINS/CNS/GPS integrated navigation based Relative filtering time  The computational times of CF and FKF are shown in Table 3. It can be seen that the time consumptions of CF and FKF are 6 3.4795 10 − × s and 5 5.5088 10 − × s, respectively. So, the timing performance of CF is better than that of FKF. There, the time consumption refers to the average time to process a frame of experimental data.

Conclusions
In this paper, a novel information fusion method based on complementary filter is proposed for the SINS/CNS/GPS integrated navigation system of an aerospace plane. The transformation algorithm of the CNS quaternion from the i frame to the n frame is designed, and the CF-based attitude/speed/position estimations are designed and implemented for SINS/CNS/GPS integrated The computational times of CF and FKF are shown in Table 3. It can be seen that the time consumptions of CF and FKF are 3.4795 × 10 −6 s and 5.5088 × 10 −5 s, respectively. So, the timing performance of CF is better than that of FKF. There, the time consumption refers to the average time to process a frame of experimental data.

Conclusions
In this paper, a novel information fusion method based on complementary filter is proposed for the SINS/CNS/GPS integrated navigation system of an aerospace plane. The transformation algorithm of the CNS quaternion from the i frame to the n frame is designed, and the CF-based attitude/speed/position estimations are designed and implemented for SINS/CNS/GPS integrated navigation system, respectively. By utilizing the time-frequency domain analysis, the optimal cut-off frequencies of complementary filters are obtained. Compared with the FKF-based method, it can improve the real-time performance and computational efficiency of information fusion of spacecraft-borne multi-sensor integrated navigation. In order to verify the proposed CF-based information fusion algorithm, a semi-physical simulation experiment platform is designed, and the algorithm validity is verified by a semi-physical simulation experiment. The experiment results show that the CF-based integration method has the advantages of fast convergence speed, acceptable noise suppression, and acceptable navigation accuracy. By comparing the two methods of CF-based and FKF-based, it can be concluded that the accuracy of the CF-based method is close to that of the most commonly used FKF-based method; however, the CF has better real-time performance and lower computational complexity than FKF. For a SINS/CNS/GPS integrated navigation system of spacecraft, the elapsed time of CF is decreased to 5.19% of FKF, while the computational efficiency of CF is increased by 94.81%. The experiment results prove that the computational efficiency of information fusion of SINS/CNS/GPS integrated navigation is improved.