Deeply Integrated GNSS/Gyro Attitude Determination System.

Attitude determination systems based on Global Navigation Satellite Systems (GNSS) work on principle of phase interferometer, using multiple receiving antennas. They rely on a good quality of carrier phase tracking, that is not the case in real dynamic environment with low signal-to-noise ratio (SNR), for example, in a ground vehicle moving through an urban area or forest. There is still a problem in providing a GNSS attitude in such common conditions. This research is focused on improving sensitivity (i.e., the capability of providing attitude at a low SNR) and the reliability of the GNSS attitude determination system. It is contrasted with the majority of publications, where precision or computational efficiency is the main goal, but sensitivity and reliability are out of their scope. In the proposed system, sensitivity improved by using two measures: (a) tracking only phase differences instead of tracking full carrier phases-this is more sensitive due to the lower dynamics of the underlying process, and (b) using deep integration with gyroscope, where all phase differences are tracked in a vector gyro-aided loop closed on user's attitude in state vector. The algorithm synthesis is given, and simulation results are presented in this article. This shows that the minimal working SNR is lowered from 27-36 dBHz (typical) down to 20 dBHz, even with a low-cost MEMS gyroscope.


Introduction
Navigation systems are intended for the determination of coordinates, velocity and spatial orientation (attitude) that is described by rotation angles between body frame and navigation frame. For attitude determination, there are widely used inertial sensors (accelerometers and gyroscopes) comprising inertial navigation systems (INS) [1,2]. In recent decades, with the growing popularity of global navigation satellite systems (GNSS), [3][4][5] are commonly used GNSS-based attitude determination systems working on the principle of radio interferometer. The advantages of inertial sensors are autonomy and absolute immunity to interferences. Their drawbacks are errors increasing with time, the need for initial alignment and the very high cost for reasonable precision. The advantages of GNSS are a relatively low and stable attitude estimation error and good cost. On the other hand, GNSS' attitude estimation error highly depends on receiving conditions, such as signal-to-noise ratio (SNR) and the presence of multipath and interference, and also on user dynamics and the distance between receiving antennas (baseline). One of the approaches to upgrading attitude determination system efficiency is developing the inertial-satellite-navigation system (ISNS), in which signals of inertial sensors and GNSS are jointly processed. The efficiency growth in ISNS is due to the fact that INS and GNSS are mutually complementary, i.e., they work on different physical principles [6]. The integrated GNSS/INS attitude determination system potentially possesses higher filtering theory [23,24] is used, and, particularly, the extended Kalman filter (EKF) is taken as the base for integration algorithm. For synthesis purposes, the original model of received signals' carrier phases is used, which consists of two components caused by the transitional and rotational movement of the antenna system. This could synthesize the algorithm for the direct tracking of inter-antenna carrier phase differences that include rotational movement only.
The proposed deeply integrated GNSS/gyro attitude determination algorithm is recommended for implementation in navigational equipment that works in harsh receiving conditions with a low SNR, multipath, shadowing and high user dynamics in urban canyons or forest environment.

Algorithm Synthesis
For 3D attitude determination, three spaced GNSS antennas (as a minimum required number) are used, and a triad of gyroscopes. Attitude estimation algorithm is synthesized by using the statistical theory of optimal nonlinear filtering [23,25]. For this purpose, we shall formulate the statistical description of: • I/Q components at the GNSS receiver's correlator output, considered as the observations containing information about attitude hidden in carrier phase differences between antennas; • Gyro measurements; • Estimated state vector's dynamic model.
The special model is used for GNSS signal's carrier phase observed at spaced antennas. In this model, carrier phase is presented as the sum of two terms: one is defined by the transitional movement of the antenna plane along the user-satellite direction, and the other is defined by the rotational movement of the antennas plane around the chosen center of the body. This representation can exclude "transitional" terms from further synthesis, and process phase differences consisting of "rotational" terms only.
For the statistical description of gyro measurements, the random walk models are used for the sensor errors: biases, scales, misalignments. The attitude will be handled in quaternion form.
The important feature of the algorithm under development is the choice of estimated state vector. In state vector, we include quaternion representing user's attitude, gyroscope biases, scale factors and axes misalignments, and unknown phase biases in radio-frequency front-ends. The stochastic dynamic model for this state vector is developed using gyro measurements as a known function.
The deep integration approach is used for algorithm synthesis. According to this approach, optimal estimation system can be considered as vector-tracking loop consisting of a vector discriminator (where all the received GNSS signals are processed) and integration filter (that processes discriminator outputs and gyroscope measurements). Discriminator equations' derivation implies the cancellation of "transitional" terms of carrier phases. Error signals at discriminator outputs are proportional to estimation errors in inter-antenna phase differences that include only "rotational" terms of carrier phases. Vector feedback signal was formed out of predicted values of inter-antenna phase differences and corrected with estimated front-end phase biases. These values are calculated from a predicted estimate of attitude quaternion inside the integration filter. The estimate of quaternion at the output of the integration filter then transformed to estimates of conventional roll, pitch and yaw angles.

Coordinate Systems and Phase Differences
Considering the Earth-Centered Earth-Fixed reference frame OXYZ ( Figure 1) as a navigation frame, it will be denoted as ece f . Moving object has a dedicated body reference frame O c X c Y c Z c with origin O c , also denoted as rpy. The locations of GNSS antenna phase centers are denoted as A j , j = 1, 3, their coordinates are considered to be known in rpy frame. Lines A 1 A 2 , A 1 A 3 , A 2 A 3 are named as baselines. Assuming that the gyro triad is located at the point c O , its axes are aligned to rpy frame.
Therefore, the gyro measures angular rate vector in rpy frame.
The direction of signal arrival from i -th navigation satellite (NS) ( = 1, i n) will be described by vector NS,i μ , specified in ecef frame where λ is the signal's wavelength. By substituting Equation (2) to (3), we can describe phase difference in a more straight-forward form keeping in mind that is a known constant vector: From (3) and (4), it is obvious that information about attitude α is contained in relative phases  Assuming that the gyro triad is located at the point O c , its axes are aligned to rpy frame. Therefore, the gyro measures angular rate vector in rpy frame.
The direction of signal arrival from i-th navigation satellite (NS) (i = 1, n) will be described by vector µ NS,i , specified in ece f frame where , j = 1, 3 from rpy frame to ece f frame is given by the equation where α = R P Y T ; R, P, Y are the roll pitch and yaw angles respectively; C ece f rpy is the coordinate transformation matrix from rpy to ece f frame.
GNSS signals observed at the points A j , j = 1, 3 are phase-shifted with respect to the signal observed at the point O c . The phase difference of each signal, relatively, O c is where λ is the signal's wavelength. By substituting Equation (2) to (3), we can describe phase difference in a more straight-forward form keeping in mind that L rpy A j O c is a known constant vector: From (3) and (4), it is obvious that information about attitude α is contained in relative phases ψ i A j (α), j = 1, 3 of signals received by antennas A j . Hence, the general principle of attitude determination using GNSS signals is based on the observation of relative phases ψ i A j , j = 1, 3, i = 1, n in GNSS receiver and estimation attitude by using functional dependency Equation (3) in some way. The extended Kalman filer is used here for this purpose.

Gyro Measurements and Attitude Integration
Angular rate measurements at the output of three-axis strapdown gyroscope can be modeled as y ω rpy ;k = I + m g;k Ω rpy;k + b g;k + n g;k (5) where Ω rpy;k = ω r,k ω p,k ω y,k T is the true angular rate vector in rpy frame; m g;k is the matrix of axes misalignments and scale factor errors; b g;k is the bias vector; n g;k is a vector of mutually independent discrete zero-mean white Gaussian noises with variances σ 2 g . The user's attitude will be defined by quaternion q ece f rpy representing a rotation of the rpy frame axes to the ecef frame axes. We can rewrite phase difference (4) as the function of this quaternion where C q ece f rpy C ece f rpy (α). The evolution of the quaternion q ece f rpy can be described by the integration of angular rate vector Ω rpy;k and earth rotation rate with the following equation where ∆ * E is the constant quaternion representing rotation of Earth over the time interval t k−1 . . . t k ; ∆ rpy;k Ω rpy (t) is the update quaternion representing the rpy frame rotation in inertial space over time interval t k−1 . . . t k ; ⊗-is the quaternion product.

Correlation Signal Processing in GNSS Receiver
All modern GNSS receivers are contain a multi-channel correlator that carries out the multiplication and accumulation of the input signal with the quadrature reference signal as follows where s i re f ,j,I t k,l , s i re f ,j,Q t k,l are the quadrature components of reference signal; y A j ;k,l is the input signal from antenna A j , amplified and down-converted in corresponding front-end; k, l are the time indexes corresponding to the double-indexed time scale t k,l shown in Figure 2.

Gyro Measurements and Attitude Integration
Angular rate measurements at the output of three-axis strapdown gyroscope can be modeled as is a vector of mutually independent discrete zero-mean white Gaussian noises with variances σ 2 g . The user's attitude will be defined by quaternion ecef rpy q representing a rotation of the rpy frame axes to the ecef frame axes. We can rewrite phase difference (4) as the function of this quaternion where * E Δ is the constant quaternion representing rotation of Earth over the time is the update quaternion representing the rpy frame rotation in inertial space over

Correlation Signal Processing in GNSS Receiver
All modern GNSS receivers are contain a multi-channel correlator that carries out the multiplication and accumulation of the input signal with the quadrature reference signal as follows  In this diagram, the second index of time moments t k,l , t k,l+1 , t k,l+2 . . . , l = 1, M, corresponds to the time step T d -the step of fast signal processing after analog-to-digital converter, and the first index of time moments t k−1 = t k−1,0 , t k = t k,0 , t k+1 = t k+1,0 . . . corresponds to slow processing with a time step T = MT d . Rapidly changing signal waveforms y A j ;k,l need to be described using "fast" samples at the moments t k,l within interval [t k , t k+1 ]. A description of slowly changing processes, such as correlator I/Q outputs, time delays, phases, Doppler shifts and orientation, is given at moments t k . We also assume that angular rate measurements from gyro are synchronized with moments t k .
When the GNSS receiver is locked, the i-th NS signal, the delay tracking error and frequency tracking error can be considered small. We suppose that the amplitudes of received i-th NS signal at different antennas are equal A i j = A i , j = 1, 3. In this case, the correlator outputs can be expressed as [6] where ϕ i 0;k is the carrier phase tracking error at the beginning of interval [t k , t k+1 ] with respect to the point O c ; ϑ i DI;k is the bit of digital information transmitted within navigation signal that can be 0 or 1; is a Doppler shift in the signal received at the point A j due to object rotation relative to the point O c (we ; ω i D,k is the Doppler frequency shift in the reference signal; n i I,j;k , n i Q,j;k are non-correlated zero-mean white Gaussian noises with variances There is an inevitable signal delay in the GNSS receiver's signal chain (front-end) from the antenna to ADC. Because of technological differences, these delays will be different for different front-ends. Therefore, we have taken into account this phase lag as an add-on δ A j to the signal phase in front-end j. Supposing the use of Code Division Multiple Access (CDMA) signals only, an add-on δ A j will be stated as equal for all navigation signals received on the antenna A j .
Quadrature signals I i j;k , Q i j,k can be expressed in a complex form . X i j;k = I i j;k + jQ i j;k where j √ −1 is an imaginary unit. For further simplification, but without loss of generality, the origin O c of rpy frame is put to the phase center of the antenna A 1 . Then, correlator outputs for i-th NS signal and the j-th antenna ( j = 1, 3) can be rewritten in the form where . n i X,j;k = n i I,j;k + j·n i Q,j;k . In [6,19], it is shown that, even with a non-coherent tracking of the i-th NS signal, there can be computed estimations of phase differences ψ i A j ;k , j = 2, 3 by using products (11) can be considered as the measurements of phase differences ψ i A j ,ψ;k , j = 2, 3 in which the "transitional" part of the carrier phase tracking error is compensated along with the term Sensors 2020, 20, 2203 7 of 23 responsible for the reference oscillator instability that is common for all three antennas. Still, there remain unknowns, δ A 2 /A 1 and δ A 3 /A 1 , caused by non-identical front-end phase biases.

Synthesis of Deeply Integrated GNSS/Gyro Attitude Determination Algorithm
Taking into account that ψ i A 2 ;k , ψ i A 3 ;k in (11) are functions of the quaternion q ece f rpy,k (6), let us define the state vector to be estimated as where S g = m g,11 m g,22 m g,33 T -the diagonal elements of the matrix m g represented as a vector; → m g = m g,12 m g,13 m g,23 T -three off-diagonal elements of the matrix m g represented as a vector; A dynamic model of quaternion q ece f rpy,k is given by (7), where the Earth rotation quaternion ∆ * E and update quaternion ∆ rpy;k are described by formulas where ω E is the Earth rotation rate; ρ k is the rotation vector representing a rotation of the rpy frame in inertial space over the time interval t k−1 . . . t k . If this interval is small, ρ k can be approximately computed by a numerical integration of angular rate Ω rpy (t) where ξ ρ;k ≈ − T 2 n g;k + n g;k−1 . For errors → m g,k and b g,k we use a stochastic model in the form of Wiener process (random walk) where ξ bg,k is a vector of three independent white Gaussian noises with zero expectations and variances σ 2 bg ; ξ→ mg,k is a vector of three independent white Gaussian noises with zero expectations and variances σ 2 mg . A similar model is used for the front-end phase biases where ξ δ,k is a vector of independent white Gaussian noises with zero expectations and variances σ 2 δ . Equations (7), (12)-(15) describe a Markov model for the state vector x k that can be represented in the general vector form as where ξ x;k = 1 2 n g;k + n g;k−1 T is a vector of 14 independent zero-mean white Gaussian noises with variances: σ 2 g (three components), σ 2 bg (three components), σ 2 mg (six components), σ 2 δ (two components); F x k−1 , y ω rpy ;k , y ω rpy ;k−1 , ξ x;k is a nonlinear function determined by the right sides of Equations (9), (12)- (15).
It is necessary to note that gyro output y ω rpy ;k is considered to be a known function of time in the context of state dynamics (16).
According to control theory, each tracking loop must contain a discriminator and filter. Complex signals (11) can be used to obtain discriminators for biased phase differences ψ i where ψ i A j ;k is the predicted (in filter) estimate of phase difference, δ A j /A 1 is the predicted estimate of front-end j phase bias (relative to front-end 1), ω i A j ,ψ;k is the predicted estimate of phase difference rate that can be computed approximately by numerical derivation of Equation (6) with substitution q ece f rpy,k →q ece f rpy,k . Hence the computation algorithm for ω i A j ,ψ;k is following: In [20], it is shown that, if the tracking error where S i discr = A i M 2 /2 is the slope of the discrimination function; n i j,k is the noise at the discriminator output.
Using the methodology described in [6], equivalent linear observations are introduced The observations (19) are combined in a vector with the following orderỹ ψ,k = For the synthesis of the integration filter, we use equivalent linear observationsỹ ψ,k , where parameters ψ i A j ;k , δ A j /A 1 are the functions of the state vector x k , i.e., S ψ Here, S ψ x k is the observations vector function, and n k = n 1 2,k , n 1 3,k . . . n n 2,k , n n 3,k T . A vector discriminator of the following form can be associated with the observations (20) u Discr,ψ;k = u 1 Discr,ψ A 2 ;k , u 1 Sensors 2020, 20, 2203 9 of 23 Now, the observations vector (20) and the state vector x k , described by the vector Equation (16), can be substituted in general state estimation equations of the extended Kalman filter [24] x where K k is a Kalman gain matrix; E k is the a posteriori matrix of estimation covariance;Ẽ k is prediction covariance matrix for the predicted value of the estimated state vectorx k ; D ξ x ,k is the covariance matrix of dynamic disturbance noise (given in Appendix A); The formulas for the derivative matrices  (24) and (25) are given in Appendix A.
Taking into account Equations (17) and (18), Equation (23) can be rewritten in the form Finally, the proposed algorithm for attitude determination with GNSS signals and the three-axis gyro is described by Equations (11) The output attitude is represented in the form of quaternion q ece f rpy,k . To obtain conventional Euler angles of roll, pitch and yaw, we have to make the following calculations where and Here, cϕ = cos(ϕ), sϕ = sin(ϕ), cϑ = cos(ϑ), sϑ = sin(ϑ); ϕ is a geodetic latitude; ϑ is a geodetic longitude.
The block diagram of proposed system is shown in Figure 3. ; ϕ is a geodetic latitude; ϑ is a geodetic longitude. The block diagram of proposed system is shown in Figure 3.   Figure 3. Block diagram of deeply integrated GNSS/Gyro attitude determination system.

Common Signal Processing Algorithms in GNSS Receiver
In Section 2.3, it is marked that the actual input for algorithm are correlator samples (8) that are generated in the GNSS receiver using reference signals s i re f ,j,I t k,l and s i re f ,j,Q t k,l . Reference signals require estimations of delays τ i A j ,k and Doppler frequency shifts ω i D Σ ,j;k , i = 1, n, j = 2, 3 in incoming GNSS signals that can be obtained with non-coherent DLL and FLL.

Frequency Lock Loop
To simplify the Doppler shift estimation algorithm, we intentionally use the same reference Doppler frequency for all three antennas. A stochastic model is used for the Doppler shift dynamic in the form where ξ i ν;k−1 is a zero-mean white Gaussian noise process with variance D ξ ν . By defining the state vector as x i ω,k = ω i D Σ ;k ν i k we can write down (30) in vector form The filter for this process can be described by the following equation where u Discr,ω;k is the output process of frequency discriminator [6] where I i j,k ω i D Σ ,k , Q i j,k ω i D Σ ,k are defined by (8). Derivatives of these components are described in [6] as The vector of loop filter coefficients K is calculated by known formulas [6] from the required bandwidth.

Delay Lock Loop
This uses a stochastic model of delay dynamic in the form where ξ i τ;k−1 is a zero-mean white Gaussian noise with the variance D ξ τ , f 0 is a signal's carrier frequency. Here, we consider ω i D Σ ;k−1 as a partially determined process, that can be obtained from FLL output: where δ ω,k−1 is an FLL tracking error. The filter for process (35) is described by Equations [6] where ω i D Σ ;k−1 is the extrapolated Doppler shift estimation for the i-th NS produced by the tracking system (32), Equation (37) is written in supposition that the ranging code delay is approximately same for all three antennas. The derivatives in Equation (37) are computed as finite differences using early and late correlators [6,26]. The loop filter coefficient K τ is calculated by a known formula [6] from the required bandwidth.

Amplitude Estimator
Signal's amplitude A i is required for the calculation of the discriminator's slope S i discr in (21). It can be shown with (9) that metric X i j,k τ i k from (38) is equal to A i M 2 with precision to noise term. Therefore, we can estimate amplitude with the following first-order filter where β is the filter coefficient.

Simulation and Results
To evaluate the performance of the proposed system, the simulation model has been developed. The model includes user dynamics, GNSS GLONASS constellation (coordinates and velocities of satellites), 3x antenna model, three-axis MEMS gyro, correlator blocks and various tracking algorithms, as well as GNSS receiver navigation algorithms. All the processes in the model, i.e., correlator and gyro outputs, are simulated at a time step of T = 1 ms. The additional accumulation of I/Q components was performed with a variable time interval of an integer number of milliseconds.

Antenna Geometry
In the GNSS attitude determination system, the error depends on baseline length and the geometry of antenna. The antenna geometry chosen for simulation is an equal-sided triangle with a side length of 1 m. Hence, the attitude errors presented below are given for baseline length 1 m. Antenna geometry and body axes (R,P,Y) alignment are shown in Figure 4a.

User Dynamics
Simulation scenario included the rotation of an antenna system around a vertical axis with an angular rate of about 360 °/s and sinusoidal tilts around the roll axis. The corresponding dynamics of roll, pitch and yaw angles are shown in Figure 5. The angular rates in RPY frame are shown in Figure 6. There was no transitional movement in the scenario.

GNSS Constellation
The attitude error also depends on the geometry of visible GNSS constellation. Only GLONASS constellation has been implemented. All signal characteristics correspond to L1OC signals. The sky view for the visible constellation is shown in Figure 4b. The simulation scenario included eight satellites with a geometric dilution of precision (GDOP) approximately equal to 2.5.

User Dynamics
Simulation scenario included the rotation of an antenna system around a vertical axis with an angular rate of about 360 • /s and sinusoidal tilts around the roll axis. The corresponding dynamics of roll, pitch and yaw angles are shown in Figure 5. The angular rates in RPY frame are shown in Figure 6. There was no transitional movement in the scenario.

User Dynamics
Simulation scenario included the rotation of an antenna system around a vertical axis with an angular rate of about 360 °/s and sinusoidal tilts around the roll axis. The corresponding dynamics of roll, pitch and yaw angles are shown in Figure 5. The angular rates in RPY frame are shown in Figure 6. There was no transitional movement in the scenario.

Refernce Oscillator Dynamics
Commonly used TCXOs as the sources of time and frequency for GNSS receivers introduce severe dynamic disturbances for tracking loops. When it comes to the calculation of receiver's sensitivity, only a proper account for the reference oscillator's phase drift allows relevant results. Therefore, the recorded phase drift in real TCXO was used. The recording technique is described in [25]. The sample rate of phase drift process is 10 ksps. The TCXO is Morion's GK-206TK. The

User Dynamics
Simulation scenario included the rotation of an antenna system around a vertical axis with an angular rate of about 360 °/s and sinusoidal tilts around the roll axis. The corresponding dynamics of roll, pitch and yaw angles are shown in Figure 5. The angular rates in RPY frame are shown in Figure 6. There was no transitional movement in the scenario.

Refernce Oscillator Dynamics
Commonly used TCXOs as the sources of time and frequency for GNSS receivers introduce severe dynamic disturbances for tracking loops. When it comes to the calculation of receiver's sensitivity, only a proper account for the reference oscillator's phase drift allows relevant results. Therefore, the recorded phase drift in real TCXO was used. The recording technique is described in [25]. The sample rate of phase drift process is 10 ksps. The TCXO is Morion's GK-206TK. The

Refernce Oscillator Dynamics
Commonly used TCXOs as the sources of time and frequency for GNSS receivers introduce severe dynamic disturbances for tracking loops. When it comes to the calculation of receiver's sensitivity, only a proper account for the reference oscillator's phase drift allows relevant results. Therefore, the recorded phase drift in real TCXO was used. The recording technique is described in [25]. The sample rate of phase drift process is 10 ksps. The TCXO is Morion's GK-206TK. The experimental Allan deviation of this TCXO is shown in Figure 7a. The sample frequency drift process over one minute intervals is shown in Figure 7b. Frequency drift is transformed to the apparent Doppler's velocity drift with (c/ f 0 ) coefficient, where c is the speed of light and f 0 is the nominal quartz frequency.

Gyroscope Parameters
A low-cost MEMS IMU MPU9250 from InvenSense is chosen as the prototype for gyro simulation. Its error parameters are shown in Table 1.

Gyroscope Parameters
A low-cost MEMS IMU MPU9250 from InvenSense is chosen as the prototype for gyro simulation. Its error parameters are shown in Table 1.

I/Q Adaptive Accumulation
The accumulation time for I/Q components in (8) determines the time interval for equivalent observations (19) in the integration filter. The lower this interval, the lower the prediction errors and higher the attitude accuracy. On the other hand, the less accumulation time for I/Q components in (8), the lower slope S i discr and the higher effect of the noise n i j,k at the discriminator output (18). Therefore, it is necessary to adjust the accumulation time (MT d ) for specific SNR to reach an optimum between noise and prediction errors. It is supposed that the accumulation time T acc is equal to the integer number of milliseconds, i.e., T acc = N acc T, where T = 1 ms, and N acc is computed by the following empirical formula N acc = round 10 0.1(50−C/N 0 ) , where C/N 0 is the estimate of average SNR in dBHz. N acc is limited in the range of 1 . . . 100.

Integration Filter Convergence and Attitude Precision
The attitude error process showing a convergence of the integration filters is presented in Figure 8.

Integration Filter Convergence and Attitude Precision
The attitude error process showing a convergence of the integration filters is presented in Figure 8. The initial attitude errors are about 2°. The front-end phase bias errors are δ  The convergence process of relative front-end phase biases δ A 2 /A 1 and δ A 3 /A 1 is shown in Figure 9. Here, δψ 2 ≡ δ A 2 /A 1 and δψ 3 ≡ δ A 3 /A 1 . One can see that these biases are estimated from initial errors 12 • and 17 • down to 0.5 • precision. 3 A A . One can see that these biases are estimated from initial errors 12° and 17° down to 0.5° precision. Figure 9. Front-end phase bias estimation errors. Figure 9. Front-end phase bias estimation errors.

Sensitivity and Attitude Precision in Range of SNR
To evaluate the attitude precision in a range of SNR and determine the lowest working SNR (sensitivity) for the proposed system, the following scenario has been applied. The system was operated for 20 s at SNR = 45-50 dBHz to reach a stationary state. Next, a slow degradation of SNR has been performed at the rate of 0.5 dBHz per second. The attitude error process for roll, pitch and yaw angles is shown in Figure 10. The SNR estimates for visible satellites are shown in Figure 11.

Sensitivity and Attitude Precision in Range of SNR
To evaluate the attitude precision in a range of SNR and determine the lowest working SNR (sensitivity) for the proposed system, the following scenario has been applied. The system was operated for 20 s at SNR = 45-50 dBHz to reach a stationary state. Next, a slow degradation of SNR has been performed at the rate of 0.5 dBHz per second. The attitude error process for roll, pitch and yaw angles is shown in Figure 10. The SNR estimates for visible satellites are shown in Figure 11. It can be seen that attitude error growth is non-linear. When SNR falls from 45-50 down to 18-23 dBHz, the attitude error increases adequately from ±3 to ±10 arcmin. Next, errors rapidly become ±23 arcmin and continue growing to ±1°. At this error level, signal tracking is totally lost. This happens at SNR 15-18 dBHz. Hence, the sensitivity of the proposed system is about 18 dBHz. It can be seen that attitude error growth is non-linear. When SNR falls from 45-50 down to 18-23 dBHz, the attitude error increases adequately from ±3 to ±10 arcmin. Next, errors rapidly become ±23 arcmin and continue growing to ±1 • . At this error level, signal tracking is totally lost. This happens at SNR 15-18 dBHz. Hence, the sensitivity of the proposed system is about 18 dBHz. It can be seen that attitude error growth is non-linear. When SNR falls from 45-50 down to 18-23 dBHz, the attitude error increases adequately from ±3 to ±10 arcmin. Next, errors rapidly become ±23 arcmin and continue growing to ±1°. At this error level, signal tracking is totally lost. This happens at SNR 15-18 dBHz. Hence, the sensitivity of the proposed system is about 18 dBHz. The rapid growth in attitude error at 55-57 s is connected with the loss of signal tracking. This can be seen from Figure 12, where the number of tracked signals during scenario is shown.
It should be noted that a loss of tracking happens at FLL/DLL, as described in Section 2.5. Therefore, the sensitivity of the whole system is limited by non-coherent tracking loops used for the tracking of the ranging code's delay and Doppler frequency (including reference oscillator's frequency drift). The deeply integrated phase-difference tracking system is apparently more sensitive. The rapid growth in attitude error at 55-57 s is connected with the loss of signal tracking. This can be seen from Figure 12, where the number of tracked signals during scenario is shown. The process of adaptation of accumulation time acc T is shown in Figure 13. There is actually no adaptation after the 34th second, when SNR is below 30 dBHz. The Kalman filter's error covariance well reflects actual attitude errors. This allows covariance It should be noted that a loss of tracking happens at FLL/DLL, as described in Section 2.5. Therefore, the sensitivity of the whole system is limited by non-coherent tracking loops used for the tracking of the ranging code's delay and Doppler frequency (including reference oscillator's frequency drift). The deeply integrated phase-difference tracking system is apparently more sensitive.
The process of adaptation of accumulation time T acc is shown in Figure 13. There is actually no adaptation after the 34th second, when SNR is below 30 dBHz.
The Kalman filter's error covariance well reflects actual attitude errors. This allows covariance estimate to be used as a measure of the resulting attitude's precision. By using the Kalman filter's error covariance, the average attitude RMS error σ RPY is obtained as a function of SNR (C/N 0 ). This dependency is shown in Figure 14. The average attitude RMS error is introduced as σ RPY = The process of adaptation of accumulation time acc T is shown in Figure 13. There is actually no adaptation after the 34th second, when SNR is below 30 dBHz. The Kalman filter's error covariance well reflects actual attitude errors. This allows covariance estimate to be used as a measure of the resulting attitude's precision. By using the Kalman filter's error covariance, the average attitude RMS error σ RPY is obtained as a function of SNR ( This dependency is shown in Figure 14. The average attitude RMS error is introduced as

Attitude Error during GNSS Outage
When GNSS signals are not available (i.e., under the bridges or in tunnels) it is still possible to estimate attitude by using a gyroscope. This is done automatically in the integration filter. In this case attitude error will be growing due to gyroscope errors. Nevertheless, in the proposed design, gyro errors are partially estimated according to observability conditions and compensated for during GNSS outage, making the error growth as low as possible.
To evaluate the attitude growth rate during GNSS outage, a Monte-Carlo simulation of 20 tests is performed. In each test system, operated normally at SNR = 35-40 dBHz for 15 seconds, GNSS outage was simulated by setting the SNR to −100 dBHz. The outage duration is 60 seconds. Randomization of gyro errors and front-end noises for each test is applied. The resulting attitude errors are shown in Figure 15. It can be seen that yaw error grows to about ±2.5°, roll and pitch errors grow by about ±1° per minute. The average attitude RMS error grows with a rate of about 1° per minute. These errors can be considered as low for the chosen type of MEMS gyro.

Attitude Error during GNSS Outage
When GNSS signals are not available (i.e., under the bridges or in tunnels) it is still possible to estimate attitude by using a gyroscope. This is done automatically in the integration filter. In this case attitude error will be growing due to gyroscope errors. Nevertheless, in the proposed design, gyro errors are partially estimated according to observability conditions and compensated for during GNSS outage, making the error growth as low as possible.
To evaluate the attitude growth rate during GNSS outage, a Monte-Carlo simulation of 20 tests is performed. In each test system, operated normally at SNR = 35-40 dBHz for 15 s, GNSS outage was simulated by setting the SNR to −100 dBHz. The outage duration is 60 s. Randomization of gyro errors and front-end noises for each test is applied. The resulting attitude errors are shown in Figure 15. It can be seen that yaw error grows to about ±2.5 • , roll and pitch errors grow by about ±1 • per minute. The average attitude RMS error grows with a rate of about 1 • per minute. These errors can be considered as low for the chosen type of MEMS gyro.
is performed. In each test system, operated normally at SNR = 35-40 dBHz for 15 seconds, GNSS outage was simulated by setting the SNR to −100 dBHz. The outage duration is 60 seconds. Randomization of gyro errors and front-end noises for each test is applied. The resulting attitude errors are shown in Figure 15. It can be seen that yaw error grows to about ±2.5°, roll and pitch errors grow by about ±1° per minute. The average attitude RMS error grows with a rate of about 1° per minute. These errors can be considered as low for the chosen type of MEMS gyro.

Sensitivity
By designing the proposed deeply integrated algorithm for a GNSS/gyro attitude determination system, we focused on improving sensitivity, i.e., providing attitude at the lowest possible SNR. As follows from simulation results, we have achieved an attitude determination error of about ±60 arcmin (1 m baseline) at the lowest possible SNR of 18 dBHz. Comparing this result with analogous result obtained from a standard algorithm, in which carrier phases from three antennas are measured by using separate PLLs, then phase differences are calculated and used for attitude determination without using a gyroscope [6,20]. In Figure 16, the analytical curve of angular root-mean-square error (RMSE) vs. SNR for standard algorithm with 1 m baselines is shown. Compared to Figure 14

Sensitivity
By designing the proposed deeply integrated algorithm for a GNSS/gyro attitude determination system, we focused on improving sensitivity, i.e., providing attitude at the lowest possible SNR. As follows from simulation results, we have achieved an attitude determination error of about ±60 arcmin (1 m baseline) at the lowest possible SNR of 18 dBHz. Comparing this result with analogous result obtained from a standard algorithm, in which carrier phases from three antennas are measured by using separate PLLs, then phase differences are calculated and used for attitude determination without using a gyroscope [6,20]. In Figure 16, the analytical curve of angular root-mean-square error (RMSE) vs. SNR for standard algorithm with 1 m baselines is shown. Compared to Figure 14 The result is achieved due to deeply integrated gyroscope/phase-difference tracking system, in which the transitional user movement and referenced oscillator's frequency drift are removed from phase difference dynamics. Most of the remained process dynamic is compensated for with gyro. The resulting sensitivity of the whole system is limited by the non-coherent tracking loops used for The result is achieved due to deeply integrated gyroscope/phase-difference tracking system, in which the transitional user movement and referenced oscillator's frequency drift are removed from phase difference dynamics. Most of the remained process dynamic is compensated for with gyro. The resulting sensitivity of the whole system is limited by the non-coherent tracking loops used for tracking the ranging code's delay and the Doppler frequency of signals coming in the reference antenna. Hence, the further sensitivity improvement can be achieved by:

•
Using all three antennas for the tracking of common Doppler frequency and code delay (+5 dB, approximately); • Using deep integration of non-coherent tracking loops with INS and ultra-stable reference oscillator; • Using pilot GNSS signals for tracking at low SNR.

Precision
The achieved (in simulation) attitude accuracy is about ±7-8 arcmin at SNR of 35-40 dBHz with 1 m baselines between antennas and using a low-cost MEMS gyro. These errors should not be related seriously because they are obtained without accounting for multipath phase errors that are the main source of error for the GNSS attitude. We also did not consider the wind-up effect and phase pattern of GNSS antennas, which are subject to specific hardware.

Ambiguity Resolution
Taking a glance at (21), one can see that the reference value of phase difference ψ i A j ;k used in discriminator (17) is not ambiguous-it is reconstructed directly from attitude and known baseline. The error signal at discriminator output is bounded within ±1 radian because it is proportional ;k + S i discr ψ i A j ;k + δ A j /A 1 ) are not ambiguous. Hence, the ambiguity resolution procedure is required only once-for initial attitude determination and, occasionally, after long GNSS outages with more than a 1-2 min duration, while in most analogous systems, this computationally heavy procedure must be performed at each step. For initial ambiguity resolution, well-known Lambda, BC-Lambda or C-Lambda methods can be used, as well as initial alignment from accelerometers and magnetic compass.

Cycle Slips
When tracking a GNSS signal with common PLL, the cycle slip may occur because of signal interruption (due to multipath or shadowing) or because of a rapid phase jerk due to user dynamics or reference quartz oscillator stress. Cycle slip means missing a number of integer phase cycles from the observed carrier phase. This is the common problem of all GNSS RTK and attitude determination systems. In the proposed system, phase difference is analogue to carrier phase. As opposed to the general GNSS receiver design, there are no separate loops tracking phase differences in this system, therefore no cycle slip in such a loop can occur. In other words, we do not track the carrier phase, so there is no place for a cycle slip to occur. The signal degradation effects mentioned above influence the system in different ways. Due to vector tracking, the signal degradation in only one channel may lead (in worst case) to a bias of discriminator (17) output. This will cause a minor increase in attitude error, but the tracking will be stabilized due to the normal work of other channels. This makes the proposed system less vulnerable to the effects that cause cycle slips in a normal GNSS receiver.

Conclusions
This article described the development of the deeply integrated algorithm of attitude determination by using GNSS signals from three spaced antennas and a three-axis gyro. The optimal nonlinear filtering theory was used. Deep integration technology (methodology) provided a sensitivity (and hence jammer immunity) improvement. In the synthesized integration algorithm, the direct estimate of inter-antenna phase difference is used in each pair of antennas. For the transitional movement tracking, non-coherent FLL and DLL are used. The block diagram of the synthesized deeply integrated GNSS/gyro attitude determination system is provided.
The dependency of attitude estimation error on SNR is analyzed using a simulation model. It is shown that, even with a low-cost MEMS gyro, the attitude RMS error of 5 arcmin is achieved at SNR = 22-23 dBHz, 16 dB lower than that of the standard algorithm with the same baseline of 1 m. The lowest SNR at which attitude determination is maintained for the proposed system is 18 dBHz, the attitude RMS error in this case is about 20 arcmin. This SNR is equivalent to a situation when the noise jammer is present with a jammer-to-signal ratio J/S = 44 dB within a 2 MHz system bandwidth. This implies that the proposed system improves jammer immunity by 10-14 dB in comparison to a standard coherent GNSS receiver with attitude function [6,24].
Author Contributions: Conceptualization and methodology, A.P.; software and validation, A.S.; formal analysis and investigation, A.S. and A.P.; writing-original draft preparation, review and editing, A.P. and A.S.; visualization, A.S.; project administration, A.P. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Acknowledgments:
Authors are grateful to their colleague Ilya V. Korogodin for his research in the field of phase-difference tracking systems and the exploration of discriminator (17) characteristics.

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

Appendix A
Here are given the formulas for