High-rate multi-GNSS attitude determination: experiments, comparisons with inertial measurement units and applications of GNSS rotational seismology to the 2011 Tohoku Mw9.0 earthquake

High-rate GNSS positioning has been widely investigated and applied in science and engineering. We extend it to high-rate attitude determination under a multi-GNSS constellation. A series of experiments of high-rate GNSS attitude determination has been conducted on a platform with three 50 Hz geodetic receivers and two high-grade inertial measurement units (IMU). The high-rate attitude solutions are computed for each of the multi-GNSS systems and the combined constellation by either using short baselines with correct ambiguity resolution or precise point positioning (PPP) and compared with the IMU measurements. In the case of a single GNSS system, the experimental results have shown that GPS is of the best accuracy, followed by GLONASS. The results with Beidou are the noisiest. The combined multi-GNSS constellation can significantly improve the high-rate attitude solutions from any single GNSS system alone, which is, in particular, most suitable for applications to any platform in slow or quasi-static motion. However, the improvement rate could depend on proper weightings of measurements from different GNSS systems in the dynamical experiments. The accuracy of baseline-based high-rate GNSS attitude solutions remains stable over time, while that of PPP-based solutions substantially degrades with time, as theoretically expected. Within a short period of time, the PPP-based high-rate yaw solutions with the combined multi-GNSS constellation are comparable in accuracy with those computed from baselines with correct ambiguity resolution in the dynamical experiments. The attitude results from either static or dynamical experiments have shown that high-rate GNSS attitude determination is sufficiently precise to measure rotatory motions. GNSS rotational seismology is applied to the 2011 Tohoku Mw9.0 earthquake, illustrating the potential of multi-GNSS to precisely detect seismic rotatory motions.


Introduction
Attitude is required for navigation, guidance and control of objects in motion such as spacecraft, aircraft, vehicles, ships and unmanned vehicles. It can be determined either by using onboard attitude sensors such as accelerometers, magnetometers and gyros (Grewal et al 2001, Gebre-Egziabher et al 2004, Crassidis et al 2007 or satellite-based and star-tracking geometrical methods (Cohen 1992, Lu 1995, Seeber 2003 or combined solutions of onboard sensors with geometrical methods. Onboard sensors of attitude measurement can be advantageously made of small size and are very precise over a certain period of time but can erroneously drift away. As a result, precise determination of attitude often introduces uncertain bias parameters and gets involved with uncertain stochastic models to describe such parameters and uncertain kinematic models (Crassidis et al 2007).
Global navigation satellite systems (GNSS) are a precise and reliable 10-state instrument or sensor of positioning, velocity, timing and attitude, as have been repeatedly proved and demonstrated experimentally since the inception of the global positioning system (Cohen 1992, Hofmann-Wellenhof et al 1992, Parkinson and Spilker 1996, Seeber 2003. In principle, GNSS use precise geometrical carrier phase observables between GNSS satellites and the receivers fixed on a platform to determine the attitude of the platform (Ellis and Creswell 1979, Evans 1986, Van Graas and Braasch 1991, Cohen 1992, Cohen et al 1993, Lu et al 1994, Cannon and Sun 1996, Lachapelle et al 1996, Park et al 2000, Ardalan and Rezvani 2015, though code observables can be very useful at the stage of integer ambiguity resolution. Thus, unlike accelerometers, magnetometers and gyros, GNSS attitude determination is free of drift biases. Due to its geometrical nature, the acc uracy of GNSS attitude estimation depends on three basic factors: noise level of carrier phases observables, residual errors of systematic types, correct ambiguity resolution for carrier phase observables and the geometrical size of a platform (Cohen 1992, Parkinson andSpilker 1996). GNSS attitude determination has been successfully applied to a variety of science and engineering problems, for example, in marine and photogrammetric applications, flying aircraft and automatic landing, aerospace applications (Cohen 1992, Lu 1995, Cannon and Sun 1996, Lachapelle et al 1996, Juang and Huang 1997, Leite and Walter 2007, Gross et al 2012. High-rate GNSS precise positioning has attracted much attention for almost two decades and found wide applications in various areas of science and engineering. Most of engineering and sports' applications are based on GNSS precise relative positioning, often aided with sensors such as accelerometers, magnetometers and gyros. In civil engineering applications, high-rate GNSS precise positioning has been applied to monitor dynamic deformation of tall structures/buildings and bridges (Lovse et al 1995, Kijewski-Correa et al 2006, Meng et al 2007, Psimoulis et al 2008, Moschas and Stiros 2011, Moschas and Stiros 2015, Im et al 2013, Kaloop and Kim 2014, Roberts and Tang 2017. When aided with low cost micro-electro-mechanical system (MEMS), highrate GNSS precise positioning can be used to provide precise positioning for aircrafts (Bischof and Schön 2017) and to continuously measure positions, velocities, and other movement information on athletes and pedestrians in sports-related applications (Waegli andSkaloud 2009, Morrison et al 2012). In geophysical applications, high-rate GNSS precise positioning has been successfully applied to measure waveforms of earthquakes (Kouba 2003, Larson et al 2003, Genrich and Bock 2006, Grapenthin and Freymueller 2011. In the case of large earthquakes, as in the case of 2011 Tohoku Mw9.0 earthquake in Japan, no reference station can be assumed to remain unmoved during the earthquake. In this case, one will have to use high-rate GNSS precise point positioning (PPP) to measure the dynamical displacements of stations. GNSS PPP has been shown, experimentally and theoretically, to reach the accuracy level of 2-4 mm in the horizontal components and sub-centimeters in the vertical comp onent within a short period of time , Shu et al 2017 and has been successfully used for structural deformation monitoring (Yigit 2016).
Although high-rate GNSS precise positioning has been widely applied in civil engineering and earth sciences, this is not the case with high-rate GNSS attitude. The major purpose of this research is to extend the work of Xu et al (2013) to GNSS attitude measurement over a short period of time for potential applications to GNSS rotational seismology and beyond. Bearing in mind the recent advance of multi-GNSS systems, this work will be further extended to the case of multi-GNSS systems. More specifically, we will focus on the three aspects in our study. The first purpose is to investigate the performance of high-rate multi-GNSS attitude determination. Since rotational seismology has been a hot topic recently, as the second purpose, we will apply high-rate GNSS attitude determination methods to the 2011 Tohoku Mw9.0 earthquake. On the other hand, when accelerometers and/or inertial measurement units (IMU) are used with GNSS to measure dynamical motion and deformation of man-made structures, they are either simply used to provide mutual confirmation of measurement results with enhanced reliability (Meng et al 2007, Moschas and Stiros 2011 or integrated with GNSS to produce the combined solutions of positions and velocities (Roberts et al 2004, Waegli and Skaloud 2009, Bock et al 2011, Morrison et al 2012. In this latter case, it is often assumed implicitly that GNSS, accelerometers and IMU would all have the same scaling. Experiments have shown that waveforms computed from these sensors can be slightly different in scaling and should be calibrated to properly integrate data from these different types of sensors (Xu et al 2013, Bischof andSchön 2017). In the case of attitude determination, attitude solutions from IMU can be biased and drift away in the long term from GNSS attitude solutions (Cohen et al 1993, Grewal et al 2001, Gebre-Egziabher et al 2004, Crassidis et al 2007, Ardalan and Rezvani 2015. Since the bias issue can be important, theoretically and practically, in integration of multi-sensors, the third purpose of this research is to further compare high-rate GNSS attitudes with those measured with IMU over a short period of time.

Methods of multi-GNSS attitude determination
GNSS attitude determination is to solve for the attitude param eters of a moving or stationary platform by using GNSS carrier phase geometrical observables of three or more antennas installed on the platform (Cohen 1992, 1996, Lu 1995. A GNSS attitude determination system is said to be dedicated if all the antennas on the platform use the identical oscillator, otherwise, such an attitude system is said to be non-dedicated (Lu et al 1994, Lu 1995, Cannon and Sun 1996. Although attitude can be represented by using any three independent angular parameters such as Euler angles, quaternion parameters with one constraint or the attitude (direction cosine) matrix with six constraints on the elements (Cohen 1992, 1996, Lu 1995, the most often used parameters in the navigation literature are the three Euler angles of roll, pitch and yaw, which are defined in the body frame with respect to the local North (N), East (E) and Down (D) directions or the NED frame (Grewal et al 2001, Cohen 1992, 1996, Lu 1995.
GNSS attitude solutions can be obtained either by using the GNSS-derived baselines or directly from raw GNSS carrier phase observables (Cohen 1992, 1996, Lu 1995. These two methods are essentially equivalent after ambiguity resolution according to the equivalence theorem of parameters in a linear or linearized (real-valued) Gauss-Markov model (Baksalary 1984). In this latter case, one will also have to simultaneously resolve the integer ambiguity unknowns of carrier phase observables. The simultaneous least squares (LS) solution of attitude and integer ambiguity unknowns has been symbolically based on the following mixed integer linear observational model (Parkinson and Spilker 1996, Teunissen 2012a, Xu 2015: where y is an n-dimensional vector of observations, whose elements are often the double differences of carrier phase observables between two satellites and two receivers in the case of attitude determination, A and B are (n × t) and (n × m) real-valued matrices of full column rank, respectively, β is a t-dimensional real-valued non-stochastic vector, i.e. β ∈ R t , and R t is defined as the t-dimensional real-valued space. z is an m-dimensional unknown integer vector, i.e. z ∈ Z m , Z m is the m-dimensional integer space. is the error vector of the observations y, which is often assumed to be of zero mean and variance-covariance matrix W −1 σ 2 , with W being a positive definite weighting matrix and σ 2 an unknown positive scalar or the (unknown) variance of unit weight. In attitude determination, β can either consist of all the position correction unknowns of the slave antennas relative to the master/main antenna or directly stand for the unknown attitude parameters (Cohen 1992, Lu 1995, Teunissen 2012b. Very often, β may also consist of other unknowns for the correction of systematic errors such as ionospheric errors, tropospheric errors, hardware delay and/or system biases. In principle, the observational equation (1) apply mathematically to both a single GNSS system and multi-GNSS systems. Nevertheless, in the case of multi-GNSS constellations, β can include some nuisance bias parameters among different GNSS systems such as intra-system biases and different time scales. Although GNSS attitude determination has often been formulated with double difference carrier phases, Teunissen (2012b) recently proposed the idea of array-aided precise point positioning (PPP) to precisely estimate the positions of the array receivers and the attitude (Henkel 2015).
In principle, the integer LS method can be directly applied to (1) to resolve the integer parameters z, as first described for GNSS precise relative positioning in geodesy by Teunissen (1995) and further significantly improved by Chang et al (2005) and Xu et al (2012). For more mathematical reports on integer LS and reduction, the reader is referred to Lenstra et al (1982), Fincke and Pohst (1985), Schnorr and Euchner (1994) and Xu (2001Xu ( , 2006Xu ( , 2012Xu ( , 2013. Nevertheless, since the multi-antenna configuration can be measured precisely a priori, a number of baseline-constrained ambiguity resolution methods have been developed under the framework of GNSS attitude determination (Lu 1995, Wang et al 2009, Teunissen 2012a. After the ambiguity resolution, one can then go ahead to solely determine the attitude of the platform. If the baselines of a GNSS attitude system are first estimated from carrier phase observables by using precise relative positioning and/or PPP, then one can directly use them to determine the attitude of the platform. Mathematically speaking, given a number of unit vectors in both the local level NED and body frames, say (u 1 , u 2 , ..., u m ) in the NED frame and (b 1 , b 2 , ..., b m ) in the body frame, we will then solve for the unknown attitude matrix R by minimizing the following objective function: subject to the orthonormality condition: (Lu 1995), where w i is a positive weighting scalar, det{R} stands for the determinant of R and I 3 is a (3 × 3) identity matrix. If the rotation matrix is first properly parameterized with, for example, three Euler angles, the conditions (2b) and (2c) are automatically satisfied, and one can then simply solve the unconstrained minimization problem (2a) for the three Euler angles. Lu (1995) also extended the minimization model (2) to fully account for the variance-covariance matrices of both u i and b i . The optimization problem (2) was posed by Wahba (1965) under a more general condition without constraining u i and b i to unit vectors. If the condition (2c) of rotation is not considered, namely, if R of (2) is only required to be orthogonal, the ordinary LS solution of R was actually solved analytically by Schönemann (1964) in his 1964 dissertation (see also Schönemann 1966), which is also known as the orthogonal procrustes problem. If the rotation matrix is desired, the orthonormal matrix of Schönemann (1966) will have to be adjusted to satisfy the condition of rotation.
Given a set of GNSS baseline solutions resolved with each of multi-GNSS constellations and the corresponding vectors in the body frame, there exist three major classes of methods to solve the weighted attitude problem (2): the R-method, the Y-method and the q-method. The R-method is to directly solve for the rotation matrix by using singular value decomposition (SVD) (Farrell and Stuelpnagel 1966), though the solution approach by Farrell and Stuelpnagel (1966) applies to problems of arbitrary dimension. In the case of 3D-rigid motion, the closed-form solution of the rotation matrix was also given by Arun et al (1987) and Horn et al (1988). If data are highly corrupted, Umeyama (1997) showed that the analytical solution with the R-method through SVD could fail and provided an improved solution. The Y-method was first proposed by Davenport (1968) and well documented in Keat (1977). The q-method, as its name stands, is to represent the rotation by using the quaternion parameters. It was first proposed by Davenport (see Keat 1977). Different variants of the q-method includes the TRIAD algorithm (Black 1964, Markley 2002, the QUEST algorithm (Shuster and Oh 1981) and the optimal and fast quaternion estimators (Markley and Mortari 2000). However, if the vectors u i and b i are not stochastically independent and/or if their elements are of different accuracy, one cannot obtain the closed-form solution of attitude but has to numerically solve for the attitude iteratively (Kanatani and Niitsuma 2012). Furthermore, if other data from IMU sensors are available, one can construct the optimal Kalman filtering of attitude (Shuster 1990, Gebre-Egziabher et al 2004, Crassidis et al 2007, Gross et al 2012. Recently, Teunissen (2012b) proposed an array-aided PPP concept for simultaneous kinematic PPP positioning of an array and its attitude. The basic idea is to form a combined observational model, which consists of two submodels: one with single difference (SD) observables of phases and codes between satellites for the array PPP and the other with double difference (DD) observables of phases and codes for the attitude determination of the array. The unknown parameters are treated separately but the correlations of the observables between the two sub-models are fully taken into account. Henkel (2015) followed the idea of Teunissen (2012b) to integrate SD observables of phases between satellites and an inertial sensor for PPP and attitude determination. Since the major purpose of Henkel (2015) was mainly for navigation applications, the tightly integrated system of PPP and attitude was required to be of low cost. As a result, only two low-cost single frequency GNSS receivers are used together with an IMU sensor to determine the attitude of the platform.

High-rate multi-GNSS experiments of attitude determination and comparison with IMU
High-rate GNSS experiments of attitude determination with multi-GNSS constellations were carried out twice on the roof of a 16-story building inside the information science campus of Wuhan University with three Trimble Net R9 receivers/ antennas on March 22 and August 22, 2014, respectively. The schematic plot of three Trimble Net R9 receivers (Ant1, Ant2 and Ant3) for the experiments on March 22, 2014 is shown in figure 1. The three GNSS antennas are firmly fixed to the shake table (grey squares) constructed for the experiments of high-rate PPP with two almost orthogonal aluminium bars, each being of a length of about 4 m, rectangular size (8.44 cm × 4.90 cm) and thickness of about 0.18-0.20 cm, as shown in thick blue lines in figure 1. The three baselines between antennas are equal to 4.60 m (Ant1 and Ant2), 3.92 m (Ant1 and Ant3) and 2.90 m (Ant2 and Ant3), respectively, and the receivers are operational at the sampling rate of 50 Hz. For the comparative purpose, we have installed an IMU-FSAS inertial measurement unit on the shake table, which is also shown in figure 1 and supposed to measure roll, pitch and yaw at the accuracy of 29 , 29 and 43 , respectively, as specified by its maker. The sampling rate of this IMU is 200 Hz. The experiments on August 22, 2014 basically followed the same setting as that on March 22, 2014, though the baselines between antennas are slightly different and the IMU has been replaced with one made by a different maker. The major configurations of GPS, Beidou, GLONASS and Galileo satellites during the experiments of March 22 and August 22 are plotted in figure 2.
In the experiments on August 22, we did not have any difficulty in tracking GNSS satellites. However, in the experiments of March 22, we did encounter some problems in locking on satellites, in particular, Beidou and Galileo satellites. More specifically, BeiDou C02 satellite could not be tracked by Ant3. A more serious problem is that after starting the first shaking experiment (at about 16:20:00, local time), Ant3 often failed to lock on and track BeiDou satellites C04, C05, C06 and C09. Although Ant2 successfully tracked two Galileo satellites, none of these two Galileo satellites were visible to Ant1 and Ant3. As a result, Galileo satellites cannot be used and will not be reported in the following experiments on March 22. The satellite numbers and the PDOP values of each GNSS system and the combined multi-GNSS constellation at Ant1 are shown in figure 3, which are reasonably good during the experiments on March 22, with 7-9 satellites for GPS, 8-10 satellites for Beidou, and 6-8 satellites for GLONASS, respectively. The PDOP value of GPS is roughly slightly larger than 2.0 in most of the time, which is better than the other two systems Beidou and GLONASS. In the static experiments on August 22, the number of satellites for each of the multi-GNSS constellations is generally smaller than that on March 22. The number of satellites and the PDOP values will not be shown here. We should note, however, that during the dynamical experiments on March 22, Beidou system failed to determine dynamical attitudes by itself, due to the loss of tracking to five Beidou satellites at Ant3, as mentioned in the above.
Although the experiments on March 22 had more problems than those on August 22, the results showed some peculiar aspects when compared with those on August 22. More specifically, most Beidou satellites lose locking at Ant3 soon after starting the first dynamical experiment such that Beidou cannot afford to determine the attitude of the platform by itself. On the other hand, GLONASS performs reasonably well when compared with GPS. Thus, we choose to report mainly the results from these experiments in this paper. Nevertheless, we believe that the results from both experiments complement each other to help the reader better understand the advantages and potential problems of multi-GNSS attitude determination. We will also include the results of the experiments on August 22 for the purpose of comparison. In what follows, unless specified, the reported results are referred to the experiments on March 22. We should also note that the multi-GNSS attitude solution is derived by combining all the carrier phase measurements from different GNSS systems.

High-rate multi-GNSS baseline-based attitude determination: static experiments
To demonstrate the highest possible accuracy and potential problems of high-rate GNSS attitude determination, either with a single or multi-GNSS constellation, we have determined the high-rate attitudes of the platform with carrier phases in relative positioning mode by correctly resolving the integer ambiguities. More precisely, we have followed the two-step procedure to determine the attitudes of the platform by first correctly resolving the integer ambiguities of carrier phase observables, computing the baselines among the three antennas and then applying the formulation (2) to compute the rotation matrix R of attitude, with all the weightings w i in (2) set to unity. Since static and kinematic GNSS solutions in precise (relative or absolute) positioning have been well known to be of different performances, we have designed the experiments to investigate both static and kinematic attitude solutions of the platform. In the case of kinematic attitude solutions, we will compare high-rate GNSS attitudes with those output from the installed IMU on the platform, since no true values of attitudes are available as a standard reference datum.
As the first part of the experiments on March 22, we have computed the high-rate GNSS attitudes of the static platform for each of the GNSS systems, namely, GPS, Beidou and GLONASS. However, the data are processed in kinematic mode by treating the motionless platform as if it were in motion. The major purpose of this static test is to demonstrate the accuracy of high-rate GNSS attitude determination with each of the multi-GNSS systems and to see how much improvement could be achieved by combining the multi-GNSS constellations. Shown in figure 4 are the high-rate GNSS attitudes determined with each of the multi-GNSS systems, namely, GPS (green lines), Beidou (red lines), GLONASS (blue lines) and their combined multi-GNSS solutions (black lines). To clearly visualize the high-rate attitude solutions with each of the GNSS systems and the combined solutions, we have added some constant shift values to the computed attitudes of yaw, pitch and roll in figure 4. As can be seen roughly from the widths and fluctuations of the lines in figure 4, the errors of GPS are the lowest, followed by GLONASS and then Beidou. We may qualitatively conclude that GPS has the best performance of high-rate attitude determination. Beidou looks like performing better than GLONASS in the component of yaw but clearly worse in the component of roll, as can be inferred from the levels of random errors and line fluctuations of figure 4. A quantitative analysis of the error statistics of these solutions from different GNSS constellations will be given later in this subsection.
The pattern of the best performances of GPS remains basically unchanged with the static results of the experiments on August 22, which are shown in figure 5. The best results of GPS from both experiments may be explained by the smallest errors of GPS measurements on one hand and likely a smaller PDOP value on the other hand. If we compare the respective results of Beidou (red lines) and GLONASS (blue lines) in figures 4 and 5, we may find that the widths of the Beidou (red) lines are consistently larger than those of the GLONASS (blue) lines but the GLONASS results show more wave-like motions, indicating that (i) the errors of Beidou are higher than GLONASS and (ii) the GLONASS attitude solutions contain some un-modelled or unknown wave-like errors in both experiments. In the case of the multi-GNSS constellation (GPS, Beidou, GLONASS and Galileo), the combined attitude solutions significantly outperform any of the attitude solutions from any single GNSS system in the experiments of March 22, as can also be seen from figure 4. However, they are only slightly better than the attitude solutions from GPS in the experiments on August 22 (compare figure 5). The static experimental attitude solutions on both March 22 and August 22 may imply that the combined multi-GNSS attitude solutions are not necessarily always much better than any single GNSS system. Instead, both experiments may actually suggest that potential improvement of a multi-GNSS constellation depend on the random errors of measurements of different GNSS systems or proper weightings of measurements from different GNSS systems, the number of satellites for each GNSS system and the corresponding PDOP values.
To be quantitatively more precise, we have used the static high-rate GNSS attitude solutions from the experiments of March 22 and August 22 to compute the standard deviations of yaw, pitch and roll angles. To see and understand the short-and long-term accuracy behaviors of high-rate attitude solutions, we have also used the static attitude solutions to compute the standard deviations of attitudes with a different length of data ranging from 1 to 50 min from each of the multi-GNSS systems and the combined constellation, respectively, which are listed in tables 1 (March 22) and 2 (August 22), respectively. Error behaviors over a short period of time can be important in GNSS seismology, since no earthquakes rupture longer than 5 min up to the present. The results in these two tables have consistently shown that the accuracy of the high-rate attitude solutions computed with precise baselines remains almost unchanged with time, even though some larger standard deviations can sometimes be clearly observed (compare the accuracy from 20 to 50 min in the GPS baselinebased attitude solutions listed in table 2).
To correctly understand the error statistics of solutions from different GNSS systems, we limit ourselves to the static attitude results up to 5 min in tables 1 and 2. As can be seen, the GPS and Beidou solutions are of the lowest and highest levels of errors, respectively, with those of GLONASS in between, which has confirmed our qualitative analysis in the above from a quick look at and simple comparison of the widths of the corresponding lines in figures 4 and 5. More precisely, the GPS attitude solutions remain the best in all the three components of yaw, pitch and roll with all the lengths of data. If we limit ourselves to the first 5 min data, GPS is, on average, better than Beidou by a factor of about 42 percent in yaw, 48 percent in pitch and 120 percent in roll for the experiments on March 22, respectively. GPS is about 32 percent better than GLONASS in yaw but roughly the same as GLONASS in both pitch and roll. In the experiments on August 22, the average improvement rates of GPS over Beidou over the first 5 min become 130 percent in yaw, 71 percent in pitch and 129 percent in roll, respectively. GLONASS performs better than Beidou by an average of 44 percent in yaw, 64 percent in pitch and 46 percent in roll, respectively.
With the increasing length of data, the static experiments of August 22 have revealed some mixed features, as can be seen from table 2. Beidou seems to perform slightly better sometimes than GPS and GLONASS in some of the components. This phenomenon of accuracy between Beidou and GLONASS may likely be due to some unknown error sources in the GLONASS measurements that cause the wave-like motions in the GLONASS attitude solutions, as analyzed qualitatively in the above and seen from figures 4 and 5. On the other hand, table 2 has shown that the standard deviations  of GPS attitude solutions are larger than those of Beidou, when the length of data is more than 30 min. These larger standard deviations may not really imply an accuracy degradation in GPS baseline-based attitude solutions. Actually, after a careful look at figure 5, we can see a clear trend in both the yaw and pitch directions of GPS baseline-based attitude solutions, which also appears in the GPS baseline solutions of August 22, as shown in figure 6. The same trend can also be observed in the multi-GNSS solutions of attitudes and baselines in figures 5 and 6 but can hardly be identified in the Beidou and GLONASS solutions, because these latter solutions are too noisy to reveal such a small trend. The observed trend of GPS attitude solutions may more probably imply that the platform is subject to some small drift over time. Although the motion of the platform is very small, it may have been clearly observed in the baseline solutions in figure 6. In addition to the error analysis, we have also computed the mean values of yaw, pitch and roll for each of the GNSS systems and the multi-GNSS constellation with the static experimental results on March 22 and August 22, which are listed in table 3. The mean values of pitch and roll with Beidou are significantly different from those with other GNSS systems on March 22. The reason remains unclear after carefully checking the computing software, the phase and code observables.  Noticing the fact that five Beidou GEO satellites are roughly distributed symmetrically around the long bar (Ant1-Ant2) of the experiment platform, we did some extra test computations with and without Beidou GEO satellites, indicating that the Beidou GEO satellites can significantly affect the mean attitudes of pitch and roll. More precisely, the mean attitudes of pitch and roll from Beidou without the five GEO satellites change to 0.910 and 2.509 degrees on March 22, and to 1.153 and 3.540 degrees on August 22, respectively. On the other hand, the mean values of pitch and roll with GLONASS are significantly different from those with other GNSS systems on August 22. The reason may be attributed to the wave-like fluctuations of baseline-based attitude solutions of the day, as explained in the above, whose maximum amplitude is much larger than the corresponding noise levels, in particular, in the roll component.

High-rate multi-GNSS baseline-based attitude determination: dynamical experiments and comparison with IMU
As the second part of the experiments, we would like to investigate the performances of high-rate GNSS attitude determinations for a platform in motion with each of the multi-GNSS systems and the combined constellation. More precisely, we have carried out four dynamical attitude experiments by pushing, pulling and/or shaking the platform. The first three dynamical experiments are controlled on purpose by mainly manipulating the platform to be subject to motion in a particular attitude direction, namely, with the order of yaw, pitch and roll, respectively. The fourth experiment is to set the platform to move freely in any of the three attitude directions. Since the true attitude values of the platform at any time epoch are unknown in dynamical experiments, we have installed a high grade IMU on the platform for the purpose of comparison. As in the case of high-rate PPP experiments reported in Xu et al (2013), there exist the problems of scaling and timing misalignment between GNSS and IMU, which have been computed after the procedure described in Xu et al (2013 In what follows, we will report the results of baseline-based high-rate attitude determination on March 22 and August 22. We should note that due to the failure of locking to five Beidou satellites at Ant3 after starting the dynamical experiments on March 22, the number of Beidou satellites was not sufficient to determine the attitudes of the platform alone. As a result, no dynamical attitude results of Beidou alone were available on March 22.
To evaluate the dynamical performances of the high-rate multi-GNSS attitude determination, we have computed the differences between the high-rate baseline-based GNSS attitude solutions and the IMU measurements during the four dynamical experiments and used them to directly compute the roots of mean squared errors for each of the multi-GNSS systems and the combined constellation. The accuracy of the combined multi-GNSS attitude solutions is listed in column 'multi-GNSS' of table 4, which includes GPS, GLONASS and Beidou. The combined attitude solutions are generally better than those from any single GNSS system. Since locking was lost on several Beidou satellites in the dynamical experiments on March 22 and since the static results from Beidou have been found to be the noisiest, we have also computed the baseline-based attitudes for the dynamical experiments with only GPS and GLONASS by leaving Beidou out, whose accuracy is listed in column 'GPS+GLONASS' in table 4. The results without Beidou obviously are slightly better than those with Beidou. On average, the multi-GNSS attitude solutions without Beidou are more accurate than those with Beidou by 0.58, 8.63 and 6.52 percent in yaw, pitch and roll, respectively. The combined solutions with only GPS and GLONASS improve those of GPS-only solutions in accuracy by 17.64, 26.35 and 6.61 percent in yaw, pitch and roll, and those of GLONASS-only solutions by 9.71, 24.12 and 35.98 percent, respectively. However, one may see that sometimes, a negligible degradation can be found in the roll component of the first dynamical experiment, which may imply that equal weights for all the phase measurements from different GNSS systems may not be appropriate and variance component estimation should probably be incorporated in the future.
In With the static results in table 2, we can also compute the average accuracy of yaw, pitch and roll for each GNSS system on August 22. Since the accuracy in the last three rows of static data (namely, 30, 40 and 50 min in the first part of precise relative positioning) may likely indicate slight motion of the platform, as explained in the end of section 3.1, it is not included in the average computation. In general, as in the case of March 22, the baseline-based static attitude solutions are much better than that of the dynamical ones. More precisely, the average Since the four independent high-rate GNSS dynamical attitude experiments on March 22 show more or less a similar error nature, we plot the high-rate baseline-based GPS and GLONASS attitude solutions during the first dynamical experiment in figure 7, together with the IMU attitude solutions. Figure 7 has shown that when compared with the IMU attitude results, the yaw components are consistently of the least errors for both GPS and GLONASS, followed by the pitch components, as numerically confirmed by the above analysis of statistics in table 4. The roll components are of the worst performances. This observation of errors between GNSS and IMU can also be readily confirmed numerically by the accuracy results listed in columns 'roll' below precise relative positioning in table 4. A careful look at the subplots of pitch and roll components in figure 7 reveals that in the first dynamical experiment on March 22 with the yaw direction set in motion, the high-rate dynamical attitude solutions in the pitch and roll directions are significantly different from the corresponding IMU measurements. The same phenomena of poor determination of the roll angles have been consistently observed in the other three dynamical experiments as well but not shown here.
There may be three possible reasons to explain the phenom enon of significant differences between the high-rate GNSS attitude solutions and the IMU measurements in the pitch direction, and in particular, the roll direction. As the first possible reason, we recall that the three Trimble Net R9 antennas (Ant1, Ant2 and Ant3) are installed on two almost orthogonal aluminium bars, as described and illustrated in the experimental design (compare figure 1). Since the bars are of a rectangular size of 8.44 cm (horizontal) and 4.90 cm (vertical) and because of the weights of the antennas, the bars could be subject to vertical flexure during the dynamical experiments to some extent, in particular, if the vertical motions are involved. As a likely result of the aluminium bars, this configuration of antennas may favor more on the yaw measurements of the platform and better match the IMU yaw outputs on the platform than the other two measurements of pitch and roll. The flexures of aluminium bars may result in extra pitch and roll values during a dynamical experiment, which should be zero in the case of static experiments. However, the IMU has been firmly installed on a rigid plate and should certainly not contain any information on pitch and roll motions due to the flexures of aluminium bars and the weights of antennas. This may explain the small values of the IMU-measured rolls and the significant errors of the GPS and GLONASS rolls when compared with the IMU measurements. Second, as is  well known, the accuracy of the vertical component in precise relative positioning is usually 2-4 times worse than that of the horizontal components. This is also the case with our experiments, as is confirmed in table 6, which lists the standard deviations of the baselines during the first 3000 s of static period in the experiments on March 22. Since yaw is affected by the horizontal components of the baselines, it should be of the best possible accuracy (Lu et al 1994). However, because both pitch and roll are directly involved with the vertical component, their accuracy will be readily affected by the poor accuracy of the vertical components of the baselines. As the last possible reason to explain the large errors of pitch, and in particular, roll components, we may mention that the baseline projected along the roll direction is shorter (compare figure 1). Thus, the corresponding accuracy of rolls can be theoretically expected to be poorer. To give the reader an impression of the combined multi-GNSS attitude determination, we show the final multi-GNSS attitude solutions of the fourth dynamical experiment on both March 22 and August 22 in figure 8, together with the IMU attitude measurements. We can also clearly see from figure 8 that the yaw component is best determined by the multi-GNSS constellation, followed by the pitch component. The roll comp onent is reconstructed with the largest uncertainty.

High-rate multi-GNSS PPP attitude determination and comparison with IMU
Although conventional PPP is generally of accuracy at the level of centimeters in the long term (of say, hours to months) (Kouba 2003, Bahadur and Nohutcu 2018, Liu et al 2018), high-rate GNSS PPP has been shown to enable to achieve the 2-4 mm accuracy of positioning over a short period of time, which is comparable with that of precise relative positioning of short baselines through exact GNSS ambiguity resolution. The reason for the high accuracy of PPP over a short period of time is that most of random errors in the long term can be treated as systematic errors within a short period of time and can be almost completely removed or substantially reduced, as experimentally demonstrated in Xu et al (2013) (see also Yigit (2016)) and systematically analyzed by Shu et al (2017). Paziewski et al (2018) recently also demonstrated a millimetre level of acc uracy of high-rate RTK and PPP with a multi-GNSS constellation. As part of the experiments in this paper, we would like to apply high-rate multi-GNSS PPP to determine attitudes over a short period of time and investigate what accuracy PPPbased high-rate attitude solutions can achieve. More precisely, we will first compute the high-rate PPP positions of the three antennas over a short period of time, turn these positions into the baselines or vectors between the antennas and then apply the optim ization model (2) to determine the PPP-based attitude solutions. We have used the procedures described in Xu et al (2013) to determine the scaling factors and timing misalignments with the waveforms of GNSS relative positioning and the two IMUs. Since precise relative positioning is more precise for very short baselines than PPP, we will directly use the scaling factors and timing misalignment in section 3.2 to correct the IMU measurements before the PPP-based high-rate attitude solutions are compared with the IMU measurements.
As in Xu et al (2013), we use the software system PANDA to compute PPP solutions of the three antennas with convergent floating ambiguities in its current version. The basic observational equations are ionosphere-free combinations of dual frequency P-code and carrier phase observables on the basis of GPS L1 and L2, Beidou B1 and B2, GLONASS L1 and L2, and Galileo E1 and E5a frequencies. The following elevation-dependent function is used to determine the weights of ionosphere-free observables: where A e is the elevation angle. The PCO/PCV models by Schmid et al (2016) are used to correct antenna phase centers. Tidal effects are corrected after the IERS conventions 2003 and the Saastamoinen model is used to correct tropospheric delays. The PANDA software system detects cycle slips of carrier phase observables by using the geometry-free combination, Melbourne-Wübbena combination, the ionospherefree combination and the loss of lock indicator. The MGEX products are provided by Wuhan University, which include precise orbit and clock products (Guo et al 2015). In the case of GPS, the IGS final precise orbit and 5s satellite clock products are directly used in the experiments on March 22 and August 22. However, due to a limited available distribution of tracking stations, we can only use the 300s satellite clock and 15 min orbit products for the other GNSS systems (http:// mgex.igs.org/IGS_MGEX_Products.php), since the experiments were actually conducted in 2014. The unknown parameters to be estimated include epoch-wise station coordinates, epoch-wise receiver clock errors and residual zenith total delays for a pre-determined length of arc (say 2 h). The cutoff angle in the following experiments is set to 10 degrees and the convergence time in our PPP processing is about 20-40 min. For more technical details on PANDA and its application to high-rate PPP solutions, the reader is referred to Liu and Ge (2003) and Xu et al (2013). Based on the high-rate PPP positions of the three antennas fixed to the motionless platform, we have computed the highrate PPP-based static attitude solutions. As a result, we can further compute the standard deviations of the three attitude angles with the static data of 1, 3, 5, 10, 20, 30, 40 and 50 min, respectively, which have been listed in table 1 for the experiments on March 22 and table 2  within a short period of time, except for the noises of measurements, almost all other errors such as multipath, troposphere and ionosphere are of systematic nature and can be cancelled out by time difference. With the increase of time, residual errors of systematic sources tend to behave randomly. As a result, the corresponding standard deviations of PPP solutions become larger with time in the short term. Nevertheless, after some certain time, the randomness of residual errors may remain stable and so are the error levels of the corresponding PPP solutions in the long term; (iii) to compare the PPP-based attitude solutions with the static data, in general, GPS has the best accuracy performance, followed by Beidou. GLONASS performs the worst in the PPP-based attitude solutions. If limited to a short period of time up to 5 min with the static data on March 22, on average, GPS can be 1.28, 0.48 and 0.75 times better than Beidou in yaw, pitch and roll, and 2.01, 0.69 and 0.61 times better than GLONASS, respectively. On average, Beidou performs slightly better than GLONASS in yaw and pitch by 32.20 and 14.60 percent, respectively but slightly worse in roll by 10.10 percent. In the case of the static experiments on August 22, although GPS performs even much better than Beidou and GLONASS, Beidou is now negligibly worse than GLONASS in yaw but much better in pitch and roll by a factor of 0.17 and 1.29 times, respectively; and finally, (iv) the PPP-based yaw solutions look quite well and stable within a short period of time up to 5 min, which likely reflects the fact that PPP positionings in the horizontal components are much better than that in the vertical component. As in the case of baseline-based attitude determination in section 3.2, we can also clearly see from figure 9 that the errors of PPP-based pitch and roll solutions are either much larger (compare the pitch and roll components in the top three panels) or at least, comparable with the IMU measurements in the case of the roll component (compare the bottom two panels). The three possible reasons for large errors in pitch and roll may still well apply here, namely, the flexure of the aluminium bars due to the dynamical motion and the weights of antennas, the large errors in the vertical comp onent of precise GNSS positioning and a shorter baseline projected along the roll direction, as explained in section 3.2.
Since there was no problem in locking on Beidou satellites on August 22, we will now further show the PPP-based attitude results of the dynamical experiments. The corresponding accuracy statistics are also listed in table 5. Table 5 shows that all the multi-GNSS systems have the best performance in determining the yaw angle. The multi-GNSS constellation basically performs better than any single GNSS system, as can be seen from the PPP-based accuracy results in table 5.
Nevertheless, if we compare the accuracy of pitch from the multi-GNSS constellation in the third and fourth dynamical experiments with that from GPS, we can again see the importance and necessity of incorporating variance component estimation in the multi-GNSS constellation, which should be investigated in the future.
In general, the GPS PPP-based attitude solutions significantly outperform those with Beidou and GLONASS in all the four dynamical experiments on August 22 (compare  table 5). More precisely, the accuracy of GPS PPP-based results is significantly better than that of Beidou by a factor of 0.34-0.77 times in yaw, with an average value of 0.53 time, a factor of 0.96-2.65 times in pitch, with an average value of 1.88 times, and a factor of 2.06-6.81 times in roll, with an average value of 3.49 times, respectively. When compared with GLONASS, the accuracy of GPS PPP-based results is better by a maximum percentage of 77.40 for the first three dynamical experiments but significantly better in the fourth experiment by a factor of 1.75, 2.46 and 1.39 times in yaw, pitch and roll, respectively. When comparing the accuracy results of Beidou with those of GLONASS, we can see that Beidou performs slightly better than GLONASS in yaw but significantly worse by a factor of 1.04 and 2.23 times in pitch and roll, respectively.
To give the reader an impression on the general performances of each of the GNSS systems and the combined multi-GNSS constellation, we show the high-rate PPP-based attitude solutions of the fourth dynamical experiment on August 22 in figure 11, together with the corresponding IMU measurements.

Dynamical rotations from the GEONET stations during the 2011 Tohoku Mw 9.0 earthquake
Rotational seismology has been a recent topic of great interest in seismology, thanks to tremendous advance in modern technology of measurement (Igel et al 2007, Lee et al 2007, Pillet and Virieux 2007, Graizer 2009a. Four international workshops have regularly been devoted to rotational seismology, since the first workshop was held at the U.S. Geological Survey, Menlo Park in 2007. For more information on rotational seismology, the reader may refer to Lee et al (2007) and the web site www.rotational-seismology.org/ managed by the international working group on rotational seismology. Although rotational motion was believed to not exist and was not actually observed more than a century ago, Reid (1910) mentioned that small rotational motions could theoretically be possible. Damage reports of earthquakes documented in the historical literature of large earthquakes may imply that such motions could be real, as shown with evidence by the rotational shifts or destructions of vertically placed structures such as chimneys and grave tombstones (Lee et al 2007, Kozák 2009), though a mechanical interpretation of such observed rotations can be possible. Recent measurement evidences with high precision gyros can be found, for example, in Igel et al (2007), Pillet and Virieux (2007) and Graizer (2009a). Since rotational motions due to large earthquakes are small (Reid 1910) and cannot be measured with conventional seismometers without attitude information, some recent measurements and theoretical research have shown that rotations due to earthquakes can reach some arc minutes (Niazi 1986, Takeo 2009) and tilts can be in the order of a few degrees (Graizer 2009b). With the modern advance of GNSS technology, it is likely reasonable at the present stage to detect such small rotations and tilts by using a proper configuration of GNSS antennas with sufficient confidence. If a configuration of GNSS antennas is sufficiently large, the accuracy of measuring a rotational motion with GNSS can be very high. For example, assuming 10 km for baselines with an acc uracy of 5 mm, we could roughly expect that this high precision GNSS configuration can detect a rotational motion of 0.2 arc second with a confidence of 95 percent. Based on multi-GNSS constellations, a short baseline of up to 50 m with an accuracy of 1-2 mm should not be a problem technically and the corresponding rigid GNSS antenna configuration would be able to detect a rotation at the best possible accuracy of about 4.1 arc seconds. As a result, multi-GNSS rotational seismology should be technologically realistic to measure a few arc seconds level of rotational motions with a large confidence in the future, not to mention that a properly designed rigid configuration of GNSS antennas with a mean baseline up to 100 m should probably not be difficult technically either.
Nevertheless, because no GNSS antennas have been purposely installed to detect seismological rotations and because GNSS antennas currently applied in earth sciences are not rigidly connected together on the surface of the Earth, a current configuration of GNSS antennas would likely detect a mixture of small rotational motions and seismic wave motions passing through the points where antennas are installed. Even being well aware of this fact of mixed phenomena, we would take the 2011 Tohoku Mw9.0 earthquake as an example and make an attempt to compute the rotational motions with the GNSS data of GEONET. At the very least, this preliminary analysis can be qualitatively constructive, since none of conventional seismometers can be used to measure rotational motions of earthquakes, unless they are rigidly fixed together.
For this rotatory analysis, we have carefully selected some GEONET stations, which are sufficiently close together to form a GNSS configuration of attitude determination with shortest possible baselines and meanwhile are as close as possible to the hypocenter of the 2011 Tohoku Mw9.0 earthquake. As a result, we select five GNSS triangles, which are shown in figure 12 Figure 13. The baseline-based dynamical rotational motions of the five GPS triangles during the 2011 Tohoku Mw9.0 earthquake. The triangles 1-5 refer to the upper-left three panels, the upper-right three panels, the middle-left three panels, the middle-right three panels, and the lower-left three panels, respectively. coseismic value of −3.96 arc seconds in the yaw direction of triangle 3 and a maximum amplitude of 11.58 arc seconds in the roll direction of triangle 5.
It is interesting to find that the negative yaw values of four triangles in the northern part of the hypocenter indicate that the 2011 Tohoku Mw9.0 earthquake triggered this northern area to rotate anti-clockwise and the yaw value in triangle 4 in the south shows that this area with triangle 4 rotates clockwise. This pattern of rotatory motions in yaw may be physically interpreted as that the Sendai area with triangles 1, 3 and 5 has been most pulled towards the hypocenter. Except for a very small negative roll value in triangle 2, all the other rotatory motions around pitch and roll are directionally consistent, implying that the Tohoku area inclines towards the Pacific and the south. We should note, however, that since all the five triangles are certainly not rigid, the reported seismic rotations in table 7 should probably be a mixture of seismic rotatory and earthquake wave motions. Nevertheless, the rotatory results demonstrate the potential of multi-GNSS to detect rotatory seismic motions, if a sufficient large configuration of GNSS antennas is properly installed, for example, on the basis of some integrated building complex.

Concluding remarks
High-rate GNSS relative positioning has been widely applied to problems in sports (Buchheit et al 2014, Hoppe et al 2018, earth sciences (Kouba 2003, Larson et al 2003, Genrich and Bock 2006, Grapenthin and Freymueller 2011 and civil engineering (Lovse et al 1995, Psimoulis et al 2008, Moschas and Stiros 2011, Im et al 2013, Kaloop and Kim 2014, Roberts and Tang 2017. We have extended high-rate GNSS positioning to high-rate GNSS attitude determination with particular attention to the current stateof-the-art multi-GNSS systems. A series of experiments of high-rate GNSS attitude determination were conducted on a platform with three 50 Hz geodetic GNSS receivers and highgrade inertial measurement units (IMU) on March 22, 2014 and August 22, 2014. Generally speaking, the experimental results on these two days have clearly shown that the highrate GPS attitude results are of the best accuracy, followed by GLONASS. The high-rate Beidou attitude results are generally the noisiest.
As the first part of experiments, we have used the doubledifference observational model of carrier phases, fixed the integer ambiguities, computed the precise baselines and then determined the high-rate multi-GNSS attitude solutions. Both the static and dynamical experiments through the comparison with the zero true values and/or the IMU measurements have confirmed that the baseline-based method is applicable to precisely determine the high-rate multi-GNSS attitudes of the platform for each GNSS system alone and its accuracy remains quite stable with time. The high-rate baseline-based GNSS attitude solutions have also shown that yaw can be best determined. But pitch and roll are much noisier, due to the well-known poor accuracy of GNSS positioning in the vertical component and our experimental settings. The combined multi-GNSS constellation can generally significantly improve the high-rate GNSS attitude solutions, in particular, in the case of static experiments with an acc uracy of 0.01-0.02 degrees for a configuration with an average baseline of about 3-5 m. Multi-GNSS attitude determination can outperform any single GNSS system to detect rotational motion for any platform in slow or quasi-static motion. The high-rate baseline-based attitude solutions in the dynamical experiments are significantly less accurate than those in the static experiments. In some cases, the experimental attitude solutions have also clearly shown that the multi-GNSS constellation may perform slightly worse than GPS in some comp onents in the dynamical experiments, which may likely be due to improper weightings of carrier phase measurements from different GNSS systems.
We have attempted to apply the PPP-based method to determine high-rate attitude solutions under the framework of the multi-GNSS systems. Within a short period of time, PPPbased high-rate attitude solutions with the combined multi-GNSS constellation of GPS, Beidou, GLONASS and Galileo are shown to enable to achieve almost the comparable level of accuracy of baseline-based attitude solutions from any single GNSS system in the dynamical experiments. The high-rate PPP-based pitch and roll solutions perform much worse than the corresponding yaw solutions, which should theoretically be consistent with the general performance of high-rate PPP positioning, since they are involved with the less accurate vertical component of GNSS positioning. The yaw solutions can further be significantly improved with the combined multi-GNSS constellation. The poor performances of the PPPbased pitch and roll solutions may again be attributed to the relatively poor accuracy of high-rate PPP positioning in the vertical direction and our experimental settings. The experimental results have clearly demonstrated that the accuracy of high-rate PPP-based attitude solutions substantially degrades with time.
The high-rate attitude solutions from both the static and dynamical experiments have shown that baseline-and PPPbased methods are capable of determining high-rate attitudes and precisely measuring rotatory motions within a short period of time. As an application example of high-rate GNSS attitude determination in rotational seismology, we analyze the GEONET data of the 2011 Tohoku Mw9.0 earthquake, even though we are well aware of the fact that the rotatory motions are a mixture of both seismic rotatory and wave motions. The rotatory results have clearly demonstrated the potential of using a multi-GNSS constellation to precisely measure seismic rotatory motions. Actually, the accuracy of the high-rate multi-GNSS static attitude solutions from our experiments can be as high as 0.01 degrees or equivalently 1.7 × 10 −4 rad, which is sufficient to precisely detect any rotatory motion induced by earthquakes larger than M5.0 (Takeo 2009). The multi-GNSS baseline solutions in the horizontal components have also been shown to readily reach the accuracy of 1 mm. If a configuration of multi-GNSS antennas is of an average 50 m baseline, the yaw accuracy is expected to reach about 4 arc-seconds, which can easily detect rotatory motions due to large earthquakes with a high confidence. Finally, we may like to note that since an attitude determination configuration contains more information than a naive PPP, the PPP-based method of attitude determination could be theoretically further improved and more work should be done along this direction. The static and dynamical experiments have clearly shown that variance component estimation should also be incorporated in multi-GNSS attitude determination, which should be further investigated. Applications in engineering structures should be highly expected in the future.