High-Accuracy Decoupling Estimation of the Systematic Coordinate Errors of an INS and Intensified High Dynamic Star Tracker Based on the Constrained Least Squares Method

Navigation accuracy is one of the key performance indicators of an inertial navigation system (INS). Requirements for an accuracy assessment of an INS in a real work environment are exceedingly urgent because of enormous differences between real work and laboratory test environments. An attitude accuracy assessment of an INS based on the intensified high dynamic star tracker (IHDST) is particularly suitable for a real complex dynamic environment. However, the coupled systematic coordinate errors of an INS and the IHDST severely decrease the attitude assessment accuracy of an INS. Given that, a high-accuracy decoupling estimation method of the above systematic coordinate errors based on the constrained least squares (CLS) method is proposed in this paper. The reference frame of the IHDST is firstly converted to be consistent with that of the INS because their reference frames are completely different. Thereafter, the decoupling estimation model of the systematic coordinate errors is established and the CLS-based optimization method is utilized to estimate errors accurately. After compensating for error, the attitude accuracy of an INS can be assessed based on IHDST accurately. Both simulated experiments and real flight experiments of aircraft are conducted, and the experimental results demonstrate that the proposed method is effective and shows excellent performance for the attitude accuracy assessment of an INS in a real work environment.


Introduction
An inertial navigation system (INS) can provide the position, velocity, and attitude knowledge of a carrier. It is an autonomous navigation system which does not rely on any external information or radiate energy to the outside, and thus it has been widely used in the military and civil fields [1][2][3][4][5]. Navigation accuracy is always one of the key performance indicators of an INS. In order to improve navigation accuracy, the error parameters of an INS should be accurately calibrated before its official use. The commonly used laboratory calibration methods are the multi-position method and its various improved methods [6][7][8]. However, when an INS works in a real dynamic environment, particularly for aircraft or other high maneuverability carriers, there always exist enormous differences between the real work and laboratory test environments. This will make the credibility of the accuracy assessment of an INS in the laboratory test environment decrease remarkably. Given that, the requirements for an accuracy assessment of an INS in a real complex dynamic environment are exceedingly urgent. So far, high-accuracy speed and position benchmarks can be provided by global positioning system (GPS) [9,10], and thus an accuracy assessment of the speed and position of an INS can be realized by comparing the output information of the INS with that of the GPS. However, a high-accuracy attitude benchmark in a real complex dynamic environment is still unavailable.
A star tracker is used to determine the attitude of a carrier by matching observation stars in the field of view (FOV) and guide stars in the star catalogue. A star tracker can reach an attitude accuracy on the arc-seconds level, and it also has the characteristic of being drift-free [11][12][13]. Given that, a star tracker can be used as a potential benchmark for an attitude accuracy assessment of an INS. However, the traditional star tracker is only suitable for approximate static conditions. Under dynamic conditions, a star spot continuously moves and forms a smeared star streak because of the long exposure time of the traditional star tracker. This will make the star spot energy disperse and the signal-to-noise ratio decrease, thus reducing the attitude accuracy and even failing to output the attitude information [14,15]. Fortunately, the dynamic performance of star trackers has been significantly improved in recent years with the development of photodetectors [16][17][18] and through progress in the dynamic performance-related algorithms of star trackers [19][20][21][22][23][24]. At present, the latest intensified high dynamic star tracker (IHDST) developed by the authors has reached an attitude accuracy on the arc-seconds level in a dynamic condition of up to 25 • /s [25,26]. The IHDST can be used as the benchmark for an attitude accuracy assessment of an INS, which is particularly suitable for a real complex dynamic environment.
However, the systematic coordinate errors of an INS and IHDST, including the installation error between the INS and IHDST as well as the misalignment error of the INS, severely decrease the accuracy of the above assessment method, and the effects of the above two errors are coupled. In an integrated navigation system of an INS and a star tracker, the measurement models of the above two are known. Therefore, the existing method tends to utilize the above measurement models to establish the equations of the Karman filter or other similar filters and then estimate the systematic coordinate errors [27]. However, in this study, an attitude accuracy assessment of an INS based on the IHDST is the subject of concern. The INS becomes the assessed object, whose measurement model is unknown. At this time, the existing filter estimation methods are no longer applicable. Given that, a high-accuracy decoupling estimation method of the above systematic coordinate errors based on the constrained least squares (CLS) method is proposed in this paper. This method only utilizes the attitude data of an INS and IHDST to accurately estimate the above two errors without knowledge of the measurement models of the INS and IHDST. Since INSs and IHDSTs have completely different reference frames, the reference frame of the IHDST should be firstly converted to be consistent with that of the INS. Thereafter, the decoupling estimation model of the above two errors is established, and the CLS-based optimization method is then utilized to estimate them accurately. After compensating for the above two errors, the attitude accuracy of the INS can be ultimately assessed by using the IHDST as the attitude benchmark. Moreover, since the decoupling estimation model established and the CLS-based optimization method utilized can estimate the installation and misalignment errors accurately, the proposed method is also suitable for improving the accuracy of an INS and star tracker integrated navigation system [27][28][29].
The remainder of this paper is organized as follows. The unified principle of reference frames of the IHDST and an INS is deduced in Section 2. Section 3 details the decoupling estimation model and the CLS-based optimization method, which are utilized to estimate installation and misalignment errors accurately. Both simulated experiments and real flight experiments of an attitude accuracy assessment of an INS based on the IHDST are conducted in Section 4, and the experimental results demonstrate the feasibility and effectiveness of the proposed method. Finally, conclusions are drawn in Section 5.

Unified Principle of Reference Frame
The reference frame of the IHDST is entirely different from that of an INS due to their different measurement principles. Therefore, the reference frame of the IHDST should be converted to be consistent with that of an INS before the attitude accuracy assessment of the INS. The relevant coordinate frames are firstly defined for convenience, and then the conversion principle of a reference frame is deduced based on them.

Coordinate Frame Definition
The coordinate frames involved in this paper are shown in Figure 1, and they are defined as follows: (1) The inertial frame O e X i Y i Z i (i-frame). The J2000.0 celestial coordinate frame, which is established at 12 terrestrial dynamical time on 1 January 2000, is selected as the inertial frame in this paper. The origin O e is set at the center of the earth; the Z i -axis is normal to the equatorial plane and points towards the north celestial pole; the X i -axis lies in the equatorial plane and points towards the vernal equinox at the establishment time; and the Y i -axis completes a right-handed orthogonal frame. The inertial frame is the reference frame of the IHDST.
(2) The earth fixed frame O e X e Y e Z e (e-frame). This frame is fixed with the earth and thus remains stationary relative to the earth. The origin O e is set at the center of the earth; the Z e -axis is normal to the equatorial plane and points towards the north celestial pole; the X e -axis lies in the equatorial plane and points towards the prime meridian at the observation time; and the Y e -axis completes a right-handed orthogonal frame.
(3) The navigation frame O n X n Y n Z n (n-frame). This frame is a local vertical frame and it is related to the local geographic latitude and longitude. The origin O n is set at the location of the carrier. The X n -axis points towards the north, the Y n -axis points towards the east, and the Z n -axis points downwards. The navigation frame is the reference frame of an INS.
(4) The IHDST frame O s X s Y s Z s (s-frame). The X s -and Y s -axes are parallel to the two vertical edges of the detector plane, respectively, the Z s -axis is along the boresight of the IHDST and points outwards, and the three axes satisfy the right-hand rule.
(5) The INS frame O g X g Y g Z g (g-frame). The three axes are in accordance with the sensitive directions of the gyro triad and the accelerometer triad of the INS, and they complete the right-handed frame.

Unified Principle of Reference Frame
The reference frame of the IHDST is entirely different from that of an INS due to their different measurement principles. Therefore, the reference frame of the IHDST should be converted to be consistent with that of an INS before the attitude accuracy assessment of the INS. The relevant coordinate frames are firstly defined for convenience, and then the conversion principle of a reference frame is deduced based on them.

Coordinate Frame Definition
The coordinate frames involved in this paper are shown in Figure 1, and they are defined as follows: (1) The inertial frame OeXiYiZi (i-frame). The J2000.0 celestial coordinate frame, which is established at 12 terrestrial dynamical time on 1 January 2000, is selected as the inertial frame in this paper. The origin Oe is set at the center of the earth; the Zi-axis is normal to the equatorial plane and points towards the north celestial pole; the Xi-axis lies in the equatorial plane and points towards the vernal equinox at the establishment time; and the Yi-axis completes a right-handed orthogonal frame. The inertial frame is the reference frame of the IHDST.
(2) The earth fixed frame OeXeYeZe (e-frame). This frame is fixed with the earth and thus remains stationary relative to the earth. The origin Oe is set at the center of the earth; the Ze-axis is normal to the equatorial plane and points towards the north celestial pole; the Xe-axis lies in the equatorial plane and points towards the prime meridian at the observation time; and the Ye-axis completes a right-handed orthogonal frame.
(3) The navigation frame OnXnYnZn (n-frame). This frame is a local vertical frame and it is related to the local geographic latitude and longitude. The origin On is set at the location of the carrier. The Xn-axis points towards the north, the Yn-axis points towards the east, and the Zn-axis points downwards. The navigation frame is the reference frame of an INS.
(4) The IHDST frame OsXsYsZs (s-frame). The Xs-and Ys-axes are parallel to the two vertical edges of the detector plane, respectively, the Zs-axis is along the boresight of the IHDST and points outwards, and the three axes satisfy the right-hand rule.
(5) The INS frame OgXgYgZg (g-frame). The three axes are in accordance with the sensitive directions of the gyro triad and the accelerometer triad of the INS, and they complete the right-handed frame.

Reference Frame Conversion
As previously described, the inertial and navigation frames are the reference frames of the IHDST and INS, respectively. Therefore, the reference frame of the IHDST should be converted from the i-frame to the n-frame before an attitude accuracy assessment of the INS. The reference frame conversion can be divided into two steps, namely, the first step from the i-frame to the e-frame and the second step from the e-frame to the n-frame.

Conversion from the i-frame to the e-frame
In theory, the conversion from the i-frame to the e-frame can be realized merely depending on the earth's rotation from the epoch of J2000.0 to the observation time t. However, besides the rotation movement of the earth, its rotation axis is also moving constantly in the i-frame due to the influence of other celestial bodies (e.g., the sun and the moon) as well as the irregularity of the earth itself [30]. Figure 2 shows the motion of the earth's axis. Firstly, the north celestial pole revolves around the north ecliptic pole clockwise with a radius of the obliquity ε due to the mutual motion of the equatorial plane and the ecliptic plane, and the rotation period is about 25,770 years. Given that, a small west movement about 50.290" of the vernal equinox is generated every year. This movement is referred to as the precession. Besides that, the north celestial pole also has a small periodic elliptical swinging mainly because of the gravitation of the sun, the moon, and other celestial bodies to the earth. This movement is referred to as the nutation, and its period is about 18.6 years. Lastly, the polar motion of the earth also exists, but this kind of motion is quite small and thus can be ignored in this paper.

Reference Frame Conversion
As previously described, the inertial and navigation frames are the reference frames of the IHDST and INS, respectively. Therefore, the reference frame of the IHDST should be converted from the i-frame to the n-frame before an attitude accuracy assessment of the INS. The reference frame conversion can be divided into two steps, namely, the first step from the i-frame to the e-frame and the second step from the e-frame to the n-frame.

Conversion from the i-frame to the e-frame
In theory, the conversion from the i-frame to the e-frame can be realized merely depending on the earth's rotation from the epoch of J2000.0 to the observation time t. However, besides the rotation movement of the earth, its rotation axis is also moving constantly in the i-frame due to the influence of other celestial bodies (e.g., the sun and the moon) as well as the irregularity of the earth itself [30]. Figure 2 shows the motion of the earth's axis. Firstly, the north celestial pole revolves around the north ecliptic pole clockwise with a radius of the obliquity ε due to the mutual motion of the equatorial plane and the ecliptic plane, and the rotation period is about 25,770 years. Given that, a small west movement about 50.290˝ of the vernal equinox is generated every year. This movement is referred to as the precession. Besides that, the north celestial pole also has a small periodic elliptical swinging mainly because of the gravitation of the sun, the moon, and other celestial bodies to the earth. This movement is referred to as the nutation, and its period is about 18.6 years. Lastly, the polar motion of the earth also exists, but this kind of motion is quite small and thus can be ignored in this paper. In summary, the earth's rotation as well as the precession and nutation should be considered and compensated for, thus accurately realizing the conversion from the i-frame at the epoch of J2000.0 to the e-frame at the observation time t. Firstly, the mean celestial coordinate frame OeXimYimZim (im-frame) is obtained when the precession is compensated to the i-frame. The celestial pole, the celestial equator, and the vernal equinox corresponding to the im-frame are referred to as the mean celestial pole, the mean celestial equator, and the mean vernal equinox, respectively. Subsequently, the instantaneous celestial coordinate frame OeXitYitZit (it-frame) at the observation time t can be obtained when the nutation is compensated to the im-frame. The it-frame is established by In summary, the earth's rotation as well as the precession and nutation should be considered and compensated for, thus accurately realizing the conversion from the i-frame at the epoch of J2000.0 to the e-frame at the observation time t. Firstly, the mean celestial coordinate frame O e X im Y im Z im (i m -frame) is obtained when the precession is compensated to the i-frame. The celestial pole, the celestial equator, and the vernal equinox corresponding to the i m -frame are referred to as the mean celestial pole, the mean celestial equator, and the mean vernal equinox, respectively. Subsequently, the instantaneous celestial coordinate frame O e X it Y it Z it (i t -frame) at the observation time t can be obtained when the nutation is compensated to the i m -frame. The i t -frame is established by using the instantaneous north celestial pole as the base point, the instantaneous celestial equator as the base circle, and the instantaneous vernal equinox as the principal point, respectively. Lastly, the e-frame at the observation time t can be obtained when the earth's rotation is compensated to the i t -frame.
Let P(t), N(t), and R(t) be the precession, nutation, and earth rotation matrices, respectively, and then the preceding frame conversions can be expressed as follows [30]: (1)

1) Precession Matrix
The precession matrix in Equation (1) can be entirely determined by the three rotation matrices, and its expression is as follows: where ζ A , θ A , and z A are the precession parameters of the IAU2000 precession model [30], and R Z and R Y are the matrices rotated around the Zand Y-axes, respectively. Let δ be the rotation angle, and then the matrices rotated around the X-, Y-, and Z-axes can be expressed as follows:

2) Nutation Matrix
Similarly, the nutation matrix N(t) can also be entirely determined by the three rotation matrices, and its expression is as follows [30]: where ε, ∆ψ, and ∆ε are the obliquity, the longitude nutation, and the obliquity nutation in the IAU2000B nutation model, respectively [30,31].

3) Earth Rotation Matrix
The earth rotation matrix R(t) is entirely determined by the Greenwich hour angle of the instantaneous vernal equinox (namely, β G ), and its expression is as follows [30]: In summary, the conversion from the i-frame at the epoch of J2000.0 to the e-frame at the observation time t can be realized in three steps according to Equation (1), and the total conversion expression is as follows: where C e i (t) is the total rotation matrix from the i-frame to the e-frame, and P(t), N(t), and R(t) are determined according to Equations (2)-(5).

Conversion from the e-frame to the n-frame
As shown in Figure 1, the conversion from the e-frame to the n-frame can be realized merely depending on the longitude (λ) and latitude (ϕ) of the IHDST and INS at the observation time t. The corresponding conversion relation can be expressed as where C n e (t) is the rotation matrix from the e-frame to the n-frame, and its expression, which is entirely determined by the three rotation matrices, can be written as According to Equations (6) and (7), the total conversion from the i-frame at the epoch of J2000.0 to the n-frame at the observation time t can be expressed as Let Q s i (t) be the original attitude matrix of the IHDST with respect to the i-frame at the observation time t. Then, the conversion attitude matrix Q s n (t) of the IHDST with respect to the n-frame can be obtained according to Equation (9). At this time, the reference frame of the IHDST has been converted to be consistent with that of the INS, and the corresponding conversion can be expressed as

Decoupling Estimation of the Systematic Coordinate Errors of an INS and IHDST
As mentioned earlier, the reference frame of the IHDST has been converted from the i-frame to the n-frame based on Equation (10), and the original attitude matrix Q s i (t) of the IHDST has been accordingly transformed into the attitude matrix Q s n (t) with respect to the n-frame. Meanwhile, the reference frame of the INS is the n-frame, and thus its original attitude matrix is Q g n (t). In theory, the attitude matrix Q s n (t) of the IHDST can be directly used as the benchmark to assess the accuracy of the attitude matrix Q g n (t) of the INS. However, the actual coordinate frames of the IHDST and INS (i.e., the s-frame and the g-frame) cannot be exactly the same, and there always exists an installation error between the two. The installation error is the systematic error, and it must be compensated for, thus improving the assessment accuracy. Furthermore, according to the characteristics of the INS, its attitude error is mainly comprised of a misalignment error and an inertial instrument error. The former is the systematic error, and it should be accurately estimated and compensated for. The above two systematic coordinate errors severely decrease the accuracy of the assessment method, and their effects are coupled. Given that, a high-accuracy decoupling estimation method of the above two systematic coordinate errors based on CLS is proposed in this paper. The decoupling estimation model of the above two errors is first established, and then the CLS-based optimization method is utilized to estimate them accurately. Finally, the attitude accuracy of the INS can be assessed quantitatively after compensating for the above two errors.

Decoupling Estimation Model
Let B g s be the installation error matrix from the s-frame to the g-frame. Then, the truth matrix Q g n (t) of an INS can be derived from the attitude matrix Q s n (t) of the IHDST as follows: Since the calculation of the attitude matrix is relatively complicated, it is necessary to transform the attitude matrix into its entirely equivalent attitude quaternion, thus simplifying the related intermediate calculation. Given that, the truth matrix Q g n (t) is transformed into the truth quaternion q g n (t), the original attitude matrix Q g n (t) of the INS is transformed into the original attitude quaternion q g n (t), and the installation error matrix B g s is transformed into the installation error quaternion b g s , respectively. On this basis, Equation (11) can be converted to its quaternionic expression: where "⊗" represents quaternionic multiplication. If the measurement moment t is selected as t = t i (i = 1, 2, 3, . . . , N), the corresponding loss function L a can be expressed as In theory, the optimal estimationb g s of the installation error quaternion can be obtained when L a takes its minimum value. However, besides the installation error between the INS and the IHDST, there also exists the misalignment error of the INS, and the effects of the above two errors are mutually coupled. The optimal estimationb g s cannot be obtained any more when merely depending on the loss function L a in Equation (13). Figure 3 shows the complete transformation relations of the coordinate frames. The original attitude quaternion of the INS is not q g n (t), but q g n (t), where n and n´represent the ideal navigation frame (n-frame) and the real navigation frame (n´-frame), respectively, and the misalignment error quaternion r n n (its entirely equivalent matrix is R n n ) always exists between the two frames. estimation model of the above two errors is first established, and then the CLS-based optimization method is utilized to estimate them accurately. Finally, the attitude accuracy of the INS can be assessed quantitatively after compensating for the above two errors.

Decoupling Estimation Model
Let g s B be the installation error matrix from the s-frame to the g-frame. Then, the truth matrix  ( ) g n Q t of an INS can be derived from the attitude matrix ( ) s n Q t of the IHDST as follows: Since the calculation of the attitude matrix is relatively complicated, it is necessary to transform the attitude matrix into its entirely equivalent attitude quaternion, thus simplifying the related intermediate calculation. Given that, the truth matrix  ( )  (11) can be converted to its quaternionic expression: where " ⊗ " represents quaternionic multiplication. If the measurement moment t is selected as t=ti (i=1, 2, 3, …, N), the corresponding loss function La can be expressed as 2 2 In theory, the optimal estimation ˆg , where n and n´ represent the ideal navigation frame (n-frame) and the real navigation frame (n´-frame), respectively, and the misalignment error quaternion ′ n n r (its entirely equivalent matrix is ′ n n R ) always exists between the two frames.  After compensating for misalignment error, the attitude quaternion of the INS can be expressed as q g n (t) = r n n ⊗ q g n (t).
Sensors 2017, 17, 2285 8 of 18 By substituting Equation (14) into (13), the complete loss function L b , including the misalignment error quaternion r n n and the installation error quaternion b g s , can be expressed as where q g n (t i ) = y i0 y i1 y i2 y i3 and Equations (15)- (17) represent the decoupling estimation model of the installation and misalignment errors.

CLS-Based Optimization Method
Equations (15) and (17), which are the objective function and the equality constraints respectively, constitute a typical CLS problem. By means of equivalent transformation, the equality constraints in Equation (17) can be rewritten as The objective function L b in Equation (15) is calculated, and the result is equal to the square sum of the quaternion deviations at N measurement moments. For the moment t i , the corresponding quaternion deviation can be expressed as The above quaternion deviation ∆q i contains a total of four components, each of which consists of eight cross terms. When calculating the square sum of the 4 components of ∆q i , each component can generate 8 square terms and 28 cross-product terms, and thus 4 components totally generate 32 square terms and 112 cross-product terms. The sum Σ i,sq of all the 32 square terms can be derived as follows: When calculating the sum Σ i,cro of all the 112 cross-product terms and combining the similar terms, the expression of the ultimate result can be simplified as where and d ijk = x ij y ik , j = 0, 1, 2, 3, k = 0, 1, 2, 3.
When the above results at moment t i are extended to all the N measurement moments, the expression of the objective function L b in Equation (15) can be derived as follows: where Σ i,cro1 , Σ i,cro2 , Σ i,cro3 , and Σ i,cro4 are determined by Equations (22)- (26). After the preceding calculations, the optimal estimationsb g s andr n n can be obtained simultaneously when L b in Equation (27) takes its minimum value under the condition of Equation (18). Defining the coefficient vector as [λ 1 λ 2 ] T , the Lagrange function can be constructed as [32] l Therefore,b g s andr n n should satisfy the following Lagrange condition [32]: Equation (29) contains a total of 10 sub equations. The last two are obtained by solving the partial derivatives of a Lagrange function l relative to λ k (k = 1, 2), and they are exactly the same as Equation (18). The other eight are the partial derivatives relative to b j and r j (j = 0, 1, 2, 3), and they can be reorganized as follows: whose matrix component forms are expressed as where In theory,b g s andr n n , as well as [λ 1 λ 2 ] T , can be entirely determined by combining Equations (18), (30), and (31). However, the above equations are nonlinear, thus resulting in the solution process being quite complicated. Given that, proper equivalent transformations are conducted on Equations (30) and (31) for a calculation simplification. By multiplying the coefficient λ 2 on both sides of Equation (30) and then substituting Equation (31) into it, the result can be derived as namely, Similarly, when multiplying the coefficient λ 1 on both sides of Equation (31) and then substituting Equation (30) into it, the result can be derived as The norms of B and R are equal to 1, and thus they are both four-dimensional nonzero vectors. According to Equations (36) and (37), B and R can be considered as the eigenvectors, which belong to the eigenvalue λ 1 λ 2 of the matrices P and Q, respectively. Since the Lagrange condition described in Equation (29) is only necessary, not sufficient, all of the four eigenvectors of P and Q should be solved and then used to calculate the objective function L b in Equation (27). The optimal estimationŝ b g s andr n n can be obtained by the eigenvectors B * and R * , which make L b take its minimum value. The expressions are as follows: As previously described,b g s andr n n can be transformed into their entirely equivalent attitude matricesB g s andR n n , respectively. After compensating for the installation error, the truth matrix Q g n (t) of the INS with respect to the n-frame can be derived from the attitude matrix Q s n (t) of the IHDST, which is shown in Equation (11). Meanwhile, after compensating for the misalignment error, the measurement attitude matrix Q g n (t) of the INS can also be expressed as where Q g n (t) is the original attitude matrix of the INS. Although the attitude matrix or quaternion can clearly represent the rotation relation, they do not have any dimension, and each component does not have any independent physical meaning either. Given that, the ultimate output attitude parameters are transformed into the three-axis Euler angles in degrees ( • ), which are entirely equivalent to the attitude matrix or quaternion. In this study, the three-rotation order is defined as Z-Y-X (i.e., 3-2-1), and the corresponding Euler angles are referred to as yaw (ϕ), pitch (θ), and roll (γ), respectively. Let the truth Euler angles of an INS according to Equation (11) and the measurement Euler angles of the INS according to Equation (39) be ( γ s , θ s , ϕ s ) and (r g , θ g , ϕ g ), respectively. Then, the absolute measurement errors (AMEs) of the attitude parameters of the INS based on Euler angles can be expressed as follows:

Experiments and Discussion
Both simulated and real experiments are conducted in this section to verify the feasibility and effectiveness of the proposed method.

Simulated Experiments
In order to verify the accuracy performance of the decoupling estimation method for systematic coordinate errors proposed in this paper, simulated experiments are conducted in this section. Firstly, the truth Euler angles of the installation error (γ B , θ B , ϕ B ) and the misalignment error (γ R , θ R , ϕ R ) are randomly generated as (0.4572, −0.0146, 0.3003) and (−0.0782, 0.4157, 0.2922) in degrees ( • ), respectively. Then, the initial attitude q s n (t 0 ) of the IHDST is randomly generated, and the subsequent attitude q s n (t) of the IHDST can be derived according to the preset motion condition. The sampling frequency of the IHDST is set as 25 Hz, and the total simulated time is set as 300 s. Subsequently, According to Equations (12) and (14), the attitude q g n (t) of the INS can be determined on the basis of q s n (t), (γ B , θ B , ϕ B ) and (γ R , θ R , ϕ R ). Lastly, Gaussian noises are added to the attitude data of the IHDST and the INS, respectively, to simulate the random error. The Gaussian noise of the IHDST is set to zero mean and 5×N arc-seconds ('') standard deviation in 1σ, while the Gaussian noise of the INS is set to zero mean and 0.01×N • standard deviation in 1σ, where N is the Noise size factor, and it is set as 1, 2, 3, 4, 5, respectively. At this time, the proposed decoupling estimation method can be utilized to estimate the installation and misalignment errors. Each set of simulations is repeated five times, and the mean values of the five times are outputted as the ultimate estimation results. either. Given that, the ultimate output attitude parameters are transformed into the three-axis Euler angles in degrees (°), which are entirely equivalent to the attitude matrix or quaternion. In this study, the three-rotation order is defined as Z-Y-X (i.e., 3-2-1), and the corresponding Euler angles are referred to as yaw (φ), pitch (θ), and roll (γ), respectively.
Let the truth Euler angles of an INS according to Equation (11)

Experiments and Discussion
Both simulated and real experiments are conducted in this section to verify the feasibility and effectiveness of the proposed method.

Simulated Experiments
In order to verify the accuracy performance of the decoupling estimation method for systematic coordinate errors proposed in this paper, simulated experiments are conducted in this section. Firstly, the truth Euler angles of the installation error (γB, θB, φB) and the misalignment error  Figure 4 and Figure 5 show the absolute values of the estimation errors of the installation and misalignment errors, respectively, which are all less than 2×10 -3 ° regardless of the size of the Gaussian noise.

Real Experiments
In order to further validate the feasibility and effectiveness of the proposed method, real assessment experiments of attitude accuracy of an INS based on the IHDST, which is self-developed in this study, are conducted under an actual flight environment of an aircraft. Figure 6 shows the self-developed IHDST, and Table 1 lists its main performance specifications.

Real Experiments
In order to further validate the feasibility and effectiveness of the proposed method, real assessment experiments of attitude accuracy of an INS based on the IHDST, which is self-developed in this study, are conducted under an actual flight environment of an aircraft. Figure 6 shows the self-developed IHDST, and Table 1 lists its main performance specifications.

Real Experiments
In order to further validate the feasibility and effectiveness of the proposed method, real assessment experiments of attitude accuracy of an INS based on the IHDST, which is self-developed in this study, are conducted under an actual flight environment of an aircraft. Figure 6 shows the self-developed IHDST, and Table 1 lists its main performance specifications.     Figure 7 shows all of the experimental equipment. The selected test aircraft is shown in Figure 7a, and the X b -, Y b -, and Z b -axes of its body frame (b-frame) are along the fuselage, wing, and vertical upward directions, respectively. The assessment setup, including the INS and IHDST, is stably mounted in the test aircraft, ensuring that the s-frame and g-frame are consistent with the b-frame as much as possible, which is shown in Figure 7b Figure  7a, and the Xb-, Yb-, and Zb-axes of its body frame (b-frame) are along the fuselage, wing, and vertical upward directions, respectively. The assessment setup, including the INS and IHDST, is stably mounted in the test aircraft, ensuring that the s-frame and g-frame are consistent with the b-frame as much as possible, which is shown in Figure 7b  After the experiments, the test aircraft landed at the airport, and then the experimental attitude data of the INS and IHDST recorded in the memory was read out and analyzed. Figure 8 shows the raw three-axis Euler angles obtained by the INS and IHDST, respectively. (γs, θs, φs) are the raw data of the IHDST with respect to the i-frame, while (γg, θg, φg) are the raw data of the INS with respect to the n-frame. The correlation characteristic between the above two groups of data cannot be directly analyzed from Figure 8. Given that, the reference frame of the IHDST should be converted from the i-frame to the n-frame according to the preceding principle in this study. The converted data of the IHDST and the raw data of the INS can be plotted in the same figure, which is shown in Figure 9. After the conversion, the above two groups of data have a significant correlation characteristic in attitude variation. After the experiments, the test aircraft landed at the airport, and then the experimental attitude data of the INS and IHDST recorded in the memory was read out and analyzed. Figure 8 shows the raw three-axis Euler angles obtained by the INS and IHDST, respectively. (γ s , θ s , ϕ s ) are the raw data of the IHDST with respect to the i-frame, while (γ g , θ g , ϕ g ) are the raw data of the INS with respect to the n-frame. The correlation characteristic between the above two groups of data cannot be directly analyzed from Figure 8. Given that, the reference frame of the IHDST should be converted from the i-frame to the n-frame according to the preceding principle in this study. The converted data of the IHDST and the raw data of the INS can be plotted in the same figure, which is shown in Figure 9. After the conversion, the above two groups of data have a significant correlation characteristic in attitude variation.  As previously described, the IHDST can be used as the benchmark for an attitude accuracy assessment of the INS. The AMEs of Euler angles of the INS before and after the proposed decoupling estimation of and compensation for the systematic coordinate errors are shown in Figure  10a,b, respectively. In Figure 10a, the AMEs of roll (γ) and pitch (θ) of the INS show significant regular fluctuation characteristics. The fluctuation amplitude is relatively large, and it is about 1°. This makes it difficult to assess the attitude errors of the INS accurately. After the proposed  As previously described, the IHDST can be used as the benchmark for an attitude accuracy assessment of the INS. The AMEs of Euler angles of the INS before and after the proposed decoupling estimation of and compensation for the systematic coordinate errors are shown in Figure  10a,b, respectively. In Figure 10a, the AMEs of roll (γ) and pitch (θ) of the INS show significant regular fluctuation characteristics. The fluctuation amplitude is relatively large, and it is about 1°. This makes it difficult to assess the attitude errors of the INS accurately. After the proposed As previously described, the IHDST can be used as the benchmark for an attitude accuracy assessment of the INS. The AMEs of Euler angles of the INS before and after the proposed decoupling estimation of and compensation for the systematic coordinate errors are shown in Figure 10a,b, respectively. In Figure 10a, the AMEs of roll (γ) and pitch (θ) of the INS show significant regular fluctuation characteristics. The fluctuation amplitude is relatively large, and it is about 1 • . This makes it difficult to assess the attitude errors of the INS accurately. After the proposed decoupling estimation of and compensation for the systematic coordinate errors, the regular fluctuation has been eliminated effectively, and the AMEs of roll (γ) and pitch (θ) of the INS are reduced by about one order of magnitude, which is shown in Figure 10b. decoupling estimation of and compensation for the systematic coordinate errors, the regular fluctuation has been eliminated effectively, and the AMEs of roll (γ) and pitch (θ) of the INS are reduced by about one order of magnitude, which is shown in Figure 10b. The mean value and standard deviation of the AMEs of the INS in Figure 10 are shown in Figure 11 and Figure 12, respectively. After the estimation of and compensation for the above two errors, the mean value and standard deviation of the AMEs of the INS have been reduced significantly. For the roll (γ) and pitch (θ), the standard deviations are reduced about 30 times after error compensation, which means the errors are more concentrated. For the yaw (φ), the mean value is reduced about 15 times after error compensation, which means the error distribution is closer to the ideal zero. The corresponding optimal estimations of the installation and misalignment errors are listed in Table 2.  The mean value and standard deviation of the AMEs of the INS in Figure 10 are shown in Figures 11 and 12, respectively. After the estimation of and compensation for the above two errors, the mean value and standard deviation of the AMEs of the INS have been reduced significantly. For the roll (γ) and pitch (θ), the standard deviations are reduced about 30 times after error compensation, which means the errors are more concentrated. For the yaw (ϕ), the mean value is reduced about 15 times after error compensation, which means the error distribution is closer to the ideal zero. The corresponding optimal estimations of the installation and misalignment errors are listed in Table 2. decoupling estimation of and compensation for the systematic coordinate errors, the regular fluctuation has been eliminated effectively, and the AMEs of roll (γ) and pitch (θ) of the INS are reduced by about one order of magnitude, which is shown in Figure 10b. The mean value and standard deviation of the AMEs of the INS in Figure 10 are shown in Figure 11 and Figure 12, respectively. After the estimation of and compensation for the above two errors, the mean value and standard deviation of the AMEs of the INS have been reduced significantly. For the roll (γ) and pitch (θ), the standard deviations are reduced about 30 times after error compensation, which means the errors are more concentrated. For the yaw (φ), the mean value is reduced about 15 times after error compensation, which means the error distribution is closer to the ideal zero. The corresponding optimal estimations of the installation and misalignment errors are listed in Table 2.

Conclusions
The attitude accuracy assessment of an INS based on the IHDST is particularly suitable for a real complex dynamic environment. The coupled systematic coordinate errors severely decrease the assessment accuracy. Given that, a high-accuracy decoupling estimation method for the above systematic coordinate errors based on CLS is proposed in this paper. This method only utilizes the attitude data of the INS and the IHDST to accurately estimate the above two errors without knowledge of the measurement models of the INS and IHDST. Both simulated and real experiments were conducted to verify the feasibility and effectiveness of the proposed method. The simulated results show that the absolute values of the estimation errors of the installation and misalignment errors are all less than 2×10 −3 ° regardless of the size of the Gaussian noise. Furthermore, the real assessment experiments of the attitude accuracy of an INS are conducted under the actual flight environment of an aircraft. The IHDST is used as the benchmark for the attitude accuracy assessment of the INS. Before error compensation, the AMEs of the roll (γ) and pitch (θ) of the INS show significant regular fluctuation characteristics and the amplitude is about 1°. By contrast, after error compensation, the regular fluctuation has been eliminated effectively, and the standard deviations of the roll (γ) and pitch (θ) are reduced about 30 times, which means the errors are more concentrated. For the yaw (φ), the mean value is reduced about 15 times after error compensation, which means the error distribution is closer to the ideal zero. Moreover, the proposed method is also suitable for improving the accuracy of an INS and star tracker integrated navigation system.

Conflicts of Interest:
The authors declare no conflict of interest.

Conclusions
The attitude accuracy assessment of an INS based on the IHDST is particularly suitable for a real complex dynamic environment. The coupled systematic coordinate errors severely decrease the assessment accuracy. Given that, a high-accuracy decoupling estimation method for the above systematic coordinate errors based on CLS is proposed in this paper. This method only utilizes the attitude data of the INS and the IHDST to accurately estimate the above two errors without knowledge of the measurement models of the INS and IHDST. Both simulated and real experiments were conducted to verify the feasibility and effectiveness of the proposed method. The simulated results show that the absolute values of the estimation errors of the installation and misalignment errors are all less than 2×10 −3 • regardless of the size of the Gaussian noise. Furthermore, the real assessment experiments of the attitude accuracy of an INS are conducted under the actual flight environment of an aircraft. The IHDST is used as the benchmark for the attitude accuracy assessment of the INS. Before error compensation, the AMEs of the roll (γ) and pitch (θ) of the INS show significant regular fluctuation characteristics and the amplitude is about 1 • . By contrast, after error compensation, the regular fluctuation has been eliminated effectively, and the standard deviations of the roll (γ) and pitch (θ) are reduced about 30 times, which means the errors are more concentrated. For the yaw (ϕ), the mean value is reduced about 15 times after error compensation, which means the error distribution is closer to the ideal zero. Moreover, the proposed method is also suitable for improving the accuracy of an INS and star tracker integrated navigation system.