Semi-analytical assessment of the relative accuracy of the GNSS/INS in railway track irregularity measurements

An aided Inertial Navigation System (INS) is increasingly exploited in precise engineering surveying, such as railway track irregularity measurement, where a high relative measurement accuracy rather than absolute accuracy is emphasized. However, how to evaluate the relative measurement accuracy of the aided INS has rarely been studied. We address this problem with a semi-analytical method to analyze the relative measurement error propagation of the Global Navigation Satellite System (GNSS) and INS integrated system, specifically for the railway track irregularity measurement application. The GNSS/INS integration in this application is simplified as a linear time-invariant stochastic system driven only by white Gaussian noise, and an analytical solution for the navigation errors in the Laplace domain is obtained by analyzing the resulting steady-state Kalman filter. Then, a time series of the error is obtained through a subsequent Monte Carlo simulation based on the derived error propagation model. The proposed analysis method is then validated through data simulation and field tests. The results indicate that a 1 mm accuracy in measuring the track irregularity is achievable for the GNSS/INS integrated system. Meanwhile, the influences of the dominant inertial sensor errors on the final measurement accuracy are analyzed quantitatively and discussed comprehensively.


Introduction
The integration of a Global Navigation Satellite System (GNSS) and an Inertial Navigation System (INS) has been widely used in weapon guidance, aviation engineering, and land mobile mapping to provide accurate georeferencing (Liu et al. , 2020;El-Sheimy and Youssef , 2020). Attention is also paid to the outstanding shortterm relative measurement accuracy of the INS in inertial surveying applications (Zhang et al. , 2013;Zhu et al. , 2019;Zhang et al. , 2020). For example, an INS aided by GNSS and an odometer has successfully applied in precise engineering surveying, such as railway track irregularity measurement and road surface roughness measurement (Chen et al. , 2015;Chen et al. , 2018;Niu et al. , 2016).
The focus and accuracy requirement are very different for the above two kinds of applications. The second kind of applications pays more attention to the absolute accuracy, which is dominated by the mid-term and longterm error components, while the railway track irregularity measurement, a typical precise inertial surveying application, is more concerned with the temporal or spatial relative measurement accuracy, like the smoothness of the estimated trajectory, as illustrated in Fig. 9 in the appendix. This difference is made more explicit in the example depicted in Fig. 1. The upper panel in this figure depicts two typical samples from the first-order Gauss-Markov processes with the same covariance but different correlation times. Considering two positioning apparatuses that are corrupted by these two error processes; in this scenario, we would conclude that they have the same accuracy for georeferencing because these two stochastic processes have the same second moment, i.e., covariance matrix. In addition, it is obvious that the correlation between the values of y(t 1 ) and y(t 2 ) is higher than that between x(t 1 ) and x(t 2 ) for any two instants t 1 and t 2 , i.e., x exhibits more rapid variations in magnitude than y. Physically, we would then expect the apparatus corrupted by error process y have better accuracy in measuring railway track irregularities, as depicted in the lower panel and discussed in the appendix. This example reveals that the covariance matrix or the propagation models of the aided INS navigation errors do not contain the information describing the temporal correlation characteristics, which actually determine the relative measurement accuracy.
For the research on precise railway track irregularity measurement by an aided INS, the following two questions are often asked 1. Question 1: Is it possible to achieve 1 mm accuracy in track irregularity identification when the carrierphase differential GNSS/INS can only provide centimeter-level accuracy? 2. Question 2: What accuracy can be expected with a Track Geometry Measuring Trolley (TGMT, as shown in Fig. 10) that uses the given GNSS/INS integrated system?
These two questions are critical for the feasibility study and system design of a TGMT based on an aided INS.
The example presented in Fig. 1 intuitively illustrates why the answer to the first question is positive, while obtaining insights into these questions requires analyzing the relative measurement error propagation of the aided INS. The previous studies were almost concentrated on the absolute accuracy analysis, for example, the propagation of the covariance matrix P e of the navigation errors, as discussed in the appendix. The time history of P e portrays how the ensemble errors accumulate with time. However, this kind of analysis cannot answer the above two questions. As discussed in the appendix, it is P e in (A.7) rather than P e that characterizes the track irregularity measurement accuracy, while this information is not contained in P e . Therefore, it is impossible to evaluate the performance of a TGMT system by studying the time history of its covariance matrix.
The main difference between the present research and previous studies lies in the following: (1) The relative measurement accuracy instead of the absolute error budget of the GNSS/INS integrated system determines its performance in railway track irregularity parameter identification. Hence, we not only care about the error budget but also are interested in the correlation characteristics of the navigation errors.
(2) A semi-analytic approach is proposed to analyze the relative measurement accuracy of the GNSS/INS system, based on which the effects of the principal inertial sensor errors on the final railway track irregularity measurement accuracy are quantitatively evaluated.
The relative measurement accuracy analysis of an aided INS system was rarely studied, but with the emergence of precise engineering surveying using the inertial technique, this analysis has received more and more attention. In our previous work Zhang et al. , 2017), we pointed out the importance of relative measurement accuracy in some mobile mapping systems and studied the short-term relative measurement accuracy of the GNSS/INS system by reading the Allan variance plots of the positioning error samples of a real GNSS/INS system.
The errors in the position, velocity and attitude solution with a GNSS/INS integrated system arise from different sources, including alignment errors, inertial sensor errors, computational inaccuracy and imperfections in navigation aids, etc. Navigation error propagation is also affected by the host vehicle trajectory. It is well known that a complete determination of GNSS/INS navigation error propagation is a complex problem, and analytical solution is available only for some extremely simple cases. Maybeck (Maybeck , 1978) analyzed the navigation error of an INS aided by position data in a simplified channel . His work suggested that the analytic or at least semi-analytic performance assessment of a GNSS/INS , where x and y have the same parameter σ but different correlation times T; the lower two are defined as �α(t) = α(t) − α(t + �t) integrated system, specific for railway track irregularity measurement, is possible, as the host vehicle, i.e., the railway track trolley, moves along a simple and previously known trajectory and experiences only low dynamic motion in the surveying procedure.
The remainder of this paper is organized as follows. We first present a GNSS/INS integration algorithm specifically designed for railway track irregularity measurement and explicitly define the railway track irregularity measurement errors. Then, the system dynamic and measurement equations are simplified according to the railway track measurement condition, resulting in a time-invariant linear system. Subsequently, the relative measurement accuracy is analyzed with a semi-analytic approach. The relative error propagation analysis is then validated by a simulation. Finally, we use the semi-analytic method to quantitatively analyze the effects of the principal inertial sensor errors on the railway track irregularity measurement.

Preliminaries
For the railway track irregularity measurement application, the GNSS carrier phase observations are postprocessed in the differential mode to provide the positions accurate to the centimeter level and then the results are fused with the Inertial Measurement Unit (IMU) data in a loosely coupled architecture. The speed measurements from odometers and a Nonholonomic Constraint (NHC) are used to enhance the navigation performance. In the following, the design of a navigation Kalman Filter (KF) for loose integration of the GNSS and INS is sketched, and the track irregularity measurement errors are explicitly defined.

Used coordinate systems
We first define the coordinate systems frequently used in this work.
• the navigation frame (n-frame): a local geographic reference frame whose origin coincides with the IMU measurement center, x-axis points toward geodetic north, z-axis is down-pointing along the ellipsoid normal , and y-axis is directed east to form a righthanded frame, i.e., the North-East-Down (NED) system. • the body frame (b-frame): an IMU body frame whose axes are the same as the IMU's body axes; it is the frame in which the accelerations and angular rates generated by the strapdown accelerometers and gyroscopes are resolved. • the vehicle frame (v-frame): the host vehicle frame, whose x-axis is along the vehicle's forward direction, z-axis points downward, and y-axis is directed outward to form a right-hand system.

System models
For the GNSS/INS integrated KF design, the error state vector is defined as In this definition, δr n = [δr N δr E δr D ] T and δv n = [δv N δv E δv D ] T are the INS-derived position and velocity errors in the n-frame, respectively, and T is the three-dimensional attitude error vector, including tilt errors and the azimuth error (Benson , 1975). δb g and δb a are the errors of the gyroscope and accelerometer biases, respectively. The gyroscope and accelerometer scale factor and cross-coupling errors are not included in this error state vector, for their influence on the final navigation solution depends on the host vehicle maneuvers. The TGMT experiences only low and weak maneuvers when it moves along the rails, and therefore the scale factor and cross-coupling error of the high-grade IMU can be safely neglected. The error state vector differential equation is written as where F is the system matrix describing the system dynamics, G is the system noise distribution matrix, and w is the system noise vector. To obtain the model for an aided INS system the time derivative of each state variable is calculated. The position, velocity and attitude error differential equations are where ω n en is the angular rate vector of the n-frame with respect to the earth frame in the n-frame and δθ = [δ cos ϕ −δϕ δ sin ϕ] T , where ϕ denotes the latitude and δϕ and δ are the errors in the latitude and longitude, respectively. C n b is the b-frame to n-frame coordinate transformation matrix; f b is the specific force vector in the b-frame, where δf b is its error vector in the b-frame and f n denotes the specific force in the n-frame; ω n ie is the earth angular rotation rate vector in the n-frame; v n is the velocity vector in the n-frame; δω n ie and δω n en denote errors of ω n ie and ω n en , respectively; (1) δg n p is the local gravity error vector in the n-frame; ω n in is the angular velocity vector, ω n in = ω n ie + ω n en with the corresponding error denoted by δω n in ; and δω b ib refers to the gyro measurement errors. More details on the above three equations can be found in (Benson , 1975;Shin , 2005).
The residual inertial sensor errors are modeled as where w a , w g are three-element vectors representing the white noise component of the accelerometers and the gyro measurements, respectively. The residual biases of the gyros and accelerometer are modeled as the firstorder Gauss-Markov process where T gb and T ab are the correlation times and w gb and w ab are the corresponding driving white noise of strength 2σ 2 /T , with σ being the root mean squared value of the process.

Measurement models
For a TGMT the available external measurements include the GNSS-derived position and the speed from the odometer and NHC. The models for these observables are given below.

GNSS position update
The position of a GNSS antenna is related to the INS position by taking into account the lever arm as follows [(Teunissen and Montenbruck , 2017), p. 831]: where r G and r I are the position vectors of the GNSS antenna phase center and the IMU measurement center, respectively, which are expressed as the latitude ϕ , l ongitude and height h. l b is the lever arm from the IMU center to the GNSS antenna phase center in the b-frame, which can be accurately measured.
, a diagonal matrix that converts the delta position in meter to delta latitude, delta longitude in radians and delta height in meter, where R M and R N are, respectively, the meridian and prime vertical radii of the curvature.
Then, the measurement innovation vector comprises the difference between the GNSS-and the INS-derived positions, and the measurement model can be derived by perturbation as where r G is the estimated position of the GNSS antenna center, r n G is the position measurement provided by the GNSS receiver, n r denotes the GNSS position measurement noise in meter, modeled as Gaussian white noise with n r ∼ N (0, R r ) , which is adequate if the GNSS sampling rate is below 1 Hz (Niu et al. , 2014).

Velocity update in the vehicle frame
The host vehicle, i.e., the TGMT, is a specifically designed trolley (as depicted in Fig. 10) that differs from civilian land vehicles since it cannot change its course arbitrarily. It can keep rigid contact with the rails in both lateral and vertical directions when moving along the track and conforms to the NHC quite well even when railway track deformation exists. More details on the TGMT can be found in (Chen et al. , 2015;Chen et al. , 2018). The trolley has only an along-track speed, which can be obtained from an odometer sensor, and the velocities in both cross-track and vertical directions are zero. Therefore, the velocity measurement in the v-frame can be written as where ṽ v wheel denotes the velocity measurement vector in the v-frame; v odo represents the along-track velocity derived from the odometer output, and n v is the velocity measurement noise with n v ∼ N (0, R v ) . Note that the noise strength of the last two components of n v are different from that of the first component and should be set according to the NHC condition.
The relationship between the velocities of the trolley wheel and the IMU is expressed as where l b wheel is the lever-arm vector from the IMU measurement center to the wheel sensor in the b-frame; C v b is the b-frame to v-frame coordinate transformation matrix, which can be computed from the misalignment angles of the v-frame with respect to the b-frame; and the estimated velocity at the wheel point is denoted by v v wheel . (10) The v-frame velocity error measurement model can be expressed as

Errors of the optimal position estimates
In the GNSS/INS integrated solution, the optimal position estimate is obtained by subtracting the optimal estimate of the INS-derived position error from the INSderived position solution as where r est (t) denotes the optimal position estimate from the aided INS, with the subscript ' est' indicating the optimal estimate; r(t) is the INS-indicated position solution, which contains errors δr ; and δ r is the optimal estimate of δr from the data-fusion KF. The error in the position estimate, denoted by δr est , is defined as The INS-derived position error is defined as where r real (t) is the true position vector of the IMU. The best estimate δ r(t) also contains an estimation error, represented as Substituting (15) and (17) into (16) yields A comparison of (18) and (19) shows that the error of the estimated position of the integrated system is exactly the estimation error of the position state variable δr(t) , i.e., δr est (t) = −δ[δr(t)] . Thus, δr est (t) is used in the following error propagation modeling.

Errors of the track irregularity estimates
As introduced in the appendix, the alignment and vertical irregularities of the railway track are both defined as differential versines in (A.2), namely, track irregularity is a relative quantity. The track irregularity measurement error with the GNSS/INS system should be evaluated by the relative navigation error as defined in (A.5). Assuming the trolley moves at a constant speed, the railway track irregularity measurement error can be computed as differential navigation errors, denoted by ▽r est : where t ∈ T , ω i ∈ � ; δr est is a stochastic process that is a function defined on the product space T × , where is a fundamental sample space and T is a subset of the real line denoting a time set of interest.
Therefore, ▽r est is also a stochastic process, and contains all necessary information to completely describe the railway track irregularity measurement accuracy of a TGMT. It is clear that the objective of the performance analysis of a TGMT based on an aided INS in measuring the track irregularity is to completely characterize the relative error process ▽r est (t, ω) . In the subsequent analysis, we study the second moment, i.e., covariance matrix, of ▽r est to evaluate the measurement performance of a given TGMT. The covariance matrix of the measurement errors is defined as with

Proposed method
Equations (3) to (5) describe the INS error dynamics using a set of nonhomogeneous differential equations with time-varying coefficients. The full determination of the INS error propagation is too complicated when taking into account the real trajectories and maneuvers. In this case, we generally use a simulation approach to aid the analysis (Titterton and Weston , 2004;Groves , 2008). The situation becomes much more complex for an aided INS because the integrated navigation solution is affected not only by INS error sources but also by the external measurement noises. Fortunately, in the railway track surveying application, the nominal trajectory and motion of the host track trolley are simple and almost deterministic, which reduces the complexity and makes the analytic assessment possible. Figure 2 depicts a flowchart of the semi-analytic assessment of the relative measurement accuracy with the aided INS in measuring railway track irregularity. The procedure consists of an analytic assessment phase (the upper part of the figure) and a Monte Carlo simulation phase (the lower part). First, the system model and measurement model are simplified according to the several assumptions that hold in the specific railway surveying scenario. Then the coupling between the channels of the INS error model vanishes, and the error state dynamics is reduced to a linear time-invariant stochastic system driven only by white Gaussian noise, as enclosed by the The Monte Carlo simulation part, enclosed in the gray box, is performed as follows: create a transfer function model based on the derived δr(s) model; then, use the simulated white noise series as the input to generate the navigation error samples, i.e., the system responses; and finally, evaluate the relative measurement accuracy based on the output navigation error samples. More details on the simulation part are given in subsequent sections.

Simplification for the railway measurement case
Assumption 1 The track trolley is assumed to move uniformly in a straight line in the south-north direction.
This assumption is based on the following facts: Highspeed railway track is usually designed with a very large radius of curvature in both the horizontal and vertical profiles and with a slope gradient smaller than 1.5%. The longest chord in railway track irregularity measurements does not exceed 300 m; thus, it is reasonable to simplify such a short curve section to a straight line and to assume that the host vehicle moves uniformly in south-north direction and remains level, in which case the attitude matrix C n b = I . As defined in appendix, the horizontal railway track deformation refers to the deviation of rails from its nominal position in the lateral and vertical directions. When the trolley moves in the southnorth direction, we can evaluate the horizontal and vertical irregularity measurement accuracy by directly analyzing the east and vertical positioning errors, respectively. Thus, this assumption would simplify the following analysis. In addition, the trolley moves at a low speed, and its attitude changes slowly, which makes the dynamics-induced errors (such as the effects of the scale factor and cross-coupling) negligible.

Assumption 2
The IMU is mounted with its sensitive axes perfectly aligned with the host vehicle axes, i.e., the b-frame is coincident with the v-frame. In this case C v b becomes the identity matrix.
This assumption is reasonable because the IMU mounting angles can be estimated and compensated for with sufficient accuracy, as discussed in (Chen et al. , 2020) Assumption 3 Lever arms of the GNSS antenna and odometer are all zeros, and GNSS positions with centimeter accuracy, obtained by the postprocessed kinematic, are available all the time, as the lever arms can be measured with sufficient accuracy and the related effect can be corrected in practice.

Assumption 4
The local gravity uncertainty is assumed to be much smaller than the accelerometer measurement error and is therefore negligible.

System model simplification
The performance of the aided INS is also influenced by the real trajectories and maneuvers. We list the key information on the trajectory, inertial sensor errors, and absolute navigation error as follows: • mean position of the trajectory: latitude = 30 • , longitude = 114 • , h = 20 m. • mean absolute positioning error: • local gravity value: 9.78 m/s 2 .
We can calculate the magnitude of each term in the INS error differential equations by substituting the parameters listed above. The error terms with a magnitude smaller than 10% of the predominant term are considered insignificant terms and thus are neglected in the  . 2 Flowchart of the semi-analytical analysis of the relative measurement accuracy of the GNSS/INS integrated system following analysis. The inertial sensor error terms δω b ib and δf b are always kept. As a result, the INS error differential equations from (3) to (5) can be simplified as where δf N , δf E , δf D are the accelerometer measurement errors in the north, east, and vertical directions, respectively. δω n ib,N , δω n ib,E , δω n ib,D are the gyroscopic measurement errors in the north, east, and vertical directions, respectively. After simplification, as shown in (23), the coupling between channels vanishes, and each channel can be analyzed separately.

Measurement model simplification
A similar simplification can be carried out for measurement Eqs. (10) and (14) by taking into account the assumptions and trajectory information listed above, yielding

Measurement error propagation model
The navigation error propagation modeling of the GNSS/ INS is analyzed in the vertical and horizontal channels. In the following, the details on the vertical channel analysis are presented, while the analysis of the horizontal channel is similar and given in the appendix.
For the vertical channel analysis, in addition to the height and vertical velocity errors, the easting gyro bias δb gE and vertical accelerometer bias δb aD should be augmented into the error state vector. Additionally, since the NHC is induced as a navigation aid, the attitude error term φ E should be augmented. Thus, the error state vector and the related state dynamics model can be written as with (23) where the subscript 'D' denotes the down direction and w D (t) represents the driving white Gaussian noise of strength Q D : where δb gE is modeled as a first-order Gauss-Markov process with correlation time T gb and driven by white Gaussian noise process w gbE . The residual bias in the vertical accelerometer δb aD is modeled as a first-order Gauss-Markov process with correlation time T ab and driven by white Gaussian noise process w abD . w aD is the white Gaussian noise that corrupted the vertical accelerometer measurement. w gE is the white Gaussian noise that corrupted the gyro measurement in the east direction.
The measurements in the KF for the aided INS in the vertical channel includes the height difference between the INS and GNSS and the vertical velocity difference between the INS and the NHC. According to (25), the velocity measurement equation can be written as where n vD is the vertical component of n v in (12) and represents the measurement uncertainty of the across-track zero velocity. According to assumption 1, the host vehicle is assumed to move in the south-north direction, in which case the east velocity shall be zero. Thus, (31) can be simplified as In this case the simplification from (31) to (32) is not valid and the measurement Eq. (31) should be used. The state dynamics equations of the continuous KF are given by It is obvious from Eqs. (28) and (34) w gbE (s) + n vD (s)

Relative measurement accuracy analysis
The position error solution δr est (s) in the Laplace domain has been obtained. Theoretically, the corresponding analytic solutions in the time domain, denoted by δr est (t) , can be obtained by an inverse Laplace transform. Then, we can completely evaluate the relative measurement accuracy of the aided INS according to (A.5) and (A.7). Unfortunately, obtaining this analytic expression δr est (t) is impossible because the inverse Laplace transformation of white noise in Eqs. (41) and (B.12) does not lead to an explicit analytic expression. Therefore, the remainder analysis should be performed with a Monte Carlo simulation to generate the samples of δr est (t) , as depicted in Fig. 2. The Monte Carlo simulation relies on repeated random sampling to obtain numerical results, which is helpful in understanding the behaviors of a stochastic system that are not amenable to analysis by the usual direct mathematical methods (Brown and Hwang , 2012). Since white Gaussian noise is a stationary and Gaussian-distributed random process with a constant spectral density function over all frequencies, if it is put into a linear time-invariant system, the system output will also be stationary and Gaussian distributed. It is evident from (41), (B.12) and (20) that ▽r est will also be a stationary, Gaussian-distributed, and ergodic random process. In that case, we would be able to conduct the performance analysis of the relative measurement accuracy with sufficiently long samples δr est (t) . The analysis procedure is as follows:  (20) and its corresponding covariance matrix P ir in (21).
The information needed for generating the input white noise are listed as follows: The first observation is that the error sequence varies within 1 cm and shows an explicit spatial correlation. Based on the generated navigation error samples, we can calculate the railway track irregularity measurement error ▽r est as introduced in the above step 4. Figure 4 is a histogram of ▽r est in the east and vertical directions, which are related to the short-wave alignment and longitudinal direction, respectively, illustrating the distribution of ▽r est . They approximately follow the Gaussian distribution, as expected. The results show that 1 mm accuracy is achievable in vertical and horizontal track irregularity measurements with a navigation-grade IMU aided by a carrier-phase differential GNSS and the NHC. This result answers the first question raised in the introduction section.

Validation and discussion
In this section, we employ a full simulation approach to validate the results obtained from the proposed semianalytical error propagation model. Figure 5 depicts the validation flowchart. The previous analysis adopted a semianalytical approach where the analytical method was aided by a simulation technique in the last step. The simulation method encompasses a system simulation where the actual filter algorithm is embedded and the effects of nonlinearities in the system are kept. The simulated raw IMU data and navigation aids corrupted by different error terms are processed by a full comprehensive aided INS algorithm, i.e., the data-fusion KF, and based on its output errors; the irregularity measurement errors are analyzed. The results from these two approaches are finally compared to evaluate the feasibility of the semi-analytical method.

Aided INS raw data simulation
The aided INS data simulation software is developed at the GNSS Research Center of Wuhan University. The software consists of an IMU data simulator, navigation aiding simulator, error source simulator, and reference navigation calculator. It can generate the outputs of three orthogonal gyros and three orthogonal accelerometers, together with the navigation aid that includes position and velocity measurements via the userdefined trajectory and motion. The perfect inertial sensor output without corruption is first derived with the inverse principle of the INS. The simulated perfect IMU data are then integrated forward with time through the INS mechanization algorithm by setting the initial navigation state to generate the true navigation solutions as the navigation reference truth. The error source simulator is designed to generate several typical stochastic processes to model the inertial sensor errors and measurement noise, including white noise, random walk, and the first-order Gauss-Markov process, by using the corresponding parameters. The simulated error terms perturbing the gyros and accelerometers are added to the perfect inertial sensor outputs to yield the noise corrupted IMU data. The noise corrupted navigation aid, such as position measurement, can be achieved by adding additive noise terms to the perfect aid samples.

Simulation settings
The host vehicle trajectory and motion settings for the aided INS data simulation are consistent with the assumptions and settings in the proceeding section for error propagation modeling. The host vehicle is assumed to move in a straight line with a constant speed in the north direction and remain level. The motion details are summarized in Table 1. At the beginning, the vehicle remains stationary for ten minutes to finish the static alignment, followed by an alternating accelerating dynamic motion to make the KF converge to a steady state in a short time. Finally, it accelerates from rest at 0.2 m/s 2 for a period of 5 s, reaching a speed of 1 m/s, and then maintains uniform motion in a straight line at this velocity for 10,000 s. Only the navigation error of the uniform motion section is studied in detail in the following sections. The error terms that perturb the onboard IMU and the navigation aids are set to the same values as those for the preceding error propagation modeling analysis. The simulated IMU dataset is sampled at 200 Hz.

Simulated data processing
The  Figure 6 shows the Short wavelength Track Irregularity (STI) measurement error distribution in both alignment (or lateral) and longitudinal-level (vertical) directions. It illustrates that the measurement errors of the GNSS/INS integrated system in the simulation case closely follow a Gaussian distribution and are consistent with the those obtained from the semi-analytical approach, as shown in Fig. 4. It should be noted that the STI measurement is used here as an example, and a similar figure can also be plotted for the Long wavelength Track Irregularity (LTI) measurement. The statistical values of the measurement errors from both the semi-analytical and simulation approaches are listed in Table 2, which demonstrates that the consistence between these two methods is better than 85%. This comparison validates the feasibility of the proposed semi-analytical approach. It is notable that the measurement errors from the semi-analytical approach are slightly greater than those in the simulation. This can be explained as follows: the coupling between channels (i.e., axes) of the INS is ignored in the error propagation modeling procedure; the original weak observability of some navigation states across an axis is reduced or lost; and the loss of state observability causes the KF to lose the opportunity for state correction, lowing the final estimation accuracy slightly. In the simulation full coupling between channels has been considered, resulting in better results.

Validation in field tests
The proposed method and the related measurement accuracy analysis can also be verified through field tests, and the results were reported in our previous work (Chen et al. , 2015;Chen et al. , 2018). The TGMT based on a GNSS/INS configuration was evaluated comprehensively for the first time in November 2013 in the newly built Lanzhou-Urumqi high speed railway. The Lanzhou-Urumqi high speed passenger railway runs from Lanzhou to Urumqi in northwestern China with a designed operation speed of 250 km/h. In the experiment, the track section of about 1000 m was surveyed. The key information on the experiment and equipment used is listed below:

• GNSS/INS integrated system: a navigation grade
Positioning and Orientation System (POS) which integrates a Ring Laser Gyro (RLG) based IMU and a high-precision NovAtel GNSS OEM6 receiver. • Independent reference: a classical TGMT using high precision total station was used to provide alignment reference values accurate to 1.4 mm. A Trimble DiNi digital level was used to obtain the referenced vertical irregularity that is accurate to about 0.3 mm. • GNSS base station receiver: Trimble NetR9 GNSS reference receiver. • GNSS antenna: NovAtel GPS-702-GGL.  mm at 98% confidence level, which are consistent with the semi-analytical and simulation results. In this figure, we take the STI result as an example, more details on the field tests and results can be found in (Chen et al. , 2018). In a field test, Dead Reckoning (DR) using the GNSS/INS integrated attitude and the travel distance is suggested to exploit the NHC's potential to the fullest and ensure the consistency between the semianalytical analysis and the field tests.

Discussion
In the error propagation modeling, the semi-analytical method is used, which combines an analytical method and a Monte Carlo simulation in the relative measurement accuracy analysis. The method is proved to be feasible and credible for the performance analysis of a TGMT based on an aided INS for specific railway track irregularity measurement applications. There is an advantage of employing the semi-analytical method over the full simulation approach. Because the error state space KF model for the aided INS is simplified and reduced to a linear time-invariant stochastic system driven by the input white Gaussian noise for the specific application, it is possible to perform the accuracy analysis beforehand without a raw IMU data generator simulation and an integrated data processing. Thus, the semi-analytical method is more computationally efficient than the full simulation approach and more convenient to evaluate the effects of any changes in sensor hardware or error sources on the final measurement accuracy. It aids the system design of a TGMT based on an aided INS without actually building the system, implementing an IMU raw data simulator, and aided INS data processing software. A limitation of this approach is the analytic portion is appropriate only for the situation that the vehicle motion is extremely simple, e.g., railway track geometry surveying. It is not a general method suitable for any ground vehicle systems.

Effect of the inertial sensor errors on the measurement accuracy
With the semi-analytical approach, it is convenient to analyze the performance of a TGMT based on GNSS/INS integration and evaluate the effects of any changes in sensor hardware and error sources on the final measurement accuracy, which is very important for sensor selection in the TGMT design procedure. Here, we present a quantitative analysis of the effects of predominant inertial sensor errors, including bias instability and measurement noise of the gyroscopes and accelerometers, on the track irregularity measurement accuracy. The bias is modeled as a first-order Gauss-Markov process characterized by the correlation time and mean squared value. Figure 8 shows measurement errors of the STI and LTI in both vertical and alignment directions versus the inertial sensor bias and random noise. It shows a significant increase in measurement errors with respect to the gyroscope bias and noise, i.e., ARW, in both the STI and LTI measurement. While there are no noticeable effects of the accelerometer bias and VRW on the measurement accuracy. The determined track position and track irregularity are highly correlated with the relative measurement accuracy of the attitude angles, including roll, pitch and heading angles. Since the TGMT motion is constrained by rails, the short term attitude accuracy This figure shows that we can evaluate the system performance with the derived error propagation model and given sensor error sources. We can conclude from the above analysis that the errors in the gyroscopic measurements are the predominant error sources for the TGMT based on GNSS/INS integration. Accelerometer errors seem to have insignificant effects on the system performance and the relative measurement accuracy. However, the accelerometer bias can give rise to roll angle measurements and finally affect the cross-level measurement. Namely, gyroscopes are critical sensors that determine the performance of the railway track irregularity measurement system. Therefore, in the system design attention should be paid to the gyroscope selection. The above result has answered question 2 raised in the introduction section.

Conclusions
In this research, the temporal and spatial relative measurement accuracy and corresponding error propagation of the GNSS/INS integration specifically for the measurement of railway track irregularities are studied. A semi-analytical method is proposed to analyze the relative measurement error propagation model. In the railway surveying, the complexity of the nonlinearity of the aided INS is reduced, allowing each channel to be analyzed separately. Considering the steady-state KF operation, true errors in the optimal estimates of an aided INS are obtained with an analytical expression in the Laplace domain. Based on the derived error propagation model, the relative measurement accuracy of a railway track geometry surveying system based on the aided INS is assessed. The proposed error model is verified by a simulation, which shows the agreement between the results with the semi-analytic and simulation method is better than 85% and 93% for the vertical and horizontal channels, respectively. Therefore, this research has answered the two fundamental questions for a TGMT based on an aided INS: 1) an accuracy of 1 mm is possible in measuring the alignment and vertical track irregularities and 2) the influence of the principal inertial sensor errors on the final measurement accuracy can be quantitatively analyzed. These conclusions are valuable for the development of a TGMT and other inertial surveying applications using an aided INS and can benefit the future research. as depicted in Fig. 9. Then, the nominal offset from the rail to the chord at mileage s, denoted by v nom (s) , can be computed from the designed geometry of the track. For example, in the straightline section, the nominal versine is zero. Actual rails do not perfectly follow the designed positions and tend to drift away from the nominal smoothness. The offset from the actual rail to the chord is called the actual versine and denoted by v real (s) , which is measured by the classic geodetic instruments or a GNSS/ INS integrated system. The difference between the real and nominal versines is defined as where �v(s) can be regarded as the versine offset from the nominal value. Then, the railway track alignment irregularity parameter is computed as It is clear from the above equation that track alignment irregularity is a relative quantity. A series of Ir, i.e., the point-to-point variation of the versine deviations, portrays the track's smoothness. In the standards on precise railway track measurement, the chord length and s are defined. The chord length can be chosen as 30 m or 300 m, and the corresponding s can then be set as 5 m or 150 m, respectively. The railway track irregularity parameters related to a 30 m chord are called short wavelength irregularities, and those related to a 300 m chord are called long wavelength irregularities. The railway track irregularities are evaluated separately in both horizontal and vertical directions, with the former called alignment and the latter called longitudinal or vertical irregularity. The vertical irregularity parameters are computed using the same method depicted in Fig. 9 but plotted in the mileage-height plane. More details on the track geometric parameters can be found in our previous work (Chen et al. , 2015;Chen et al. , 2018).

Relative measurement accuracy
The navigation error process in an aided INS's solutions is stochastic. Let be a fundamental sample space and T be a subset of the real line denoting a time set of interest. Then, the navigation error process can be defined as a real-valued function of two arguments e(t, ω) , where the first argument is an element of T and the second an element of . For any fixed t ∈ T , e(t, ·) is a random variable. If we fix the second argument, for each point ω i ∈ � there is associated with a time function e(·, ω i ) , whose value at each time instant is a sample from the stochastic process.
It is common practice to describe the stochastic navigation error by the associated first two moments, i.e., the mean value function and covariance matrix. The mean value function m e (t) and covariance matrix P e (·) of the process e(·) are defined for all t ∈ T by In the preceding section, we notice that it is the error variation between the navigation error e at any time instant t and t + t that inherently determines the track irregularity measurement accuracy. Here, we denote the relative navigation error for all t ∈ T by The corresponding mean value function and covariance matrix are defined as Obviously, the information in (A.7), which characterize the TGMT's relative measurement accuracy, is not available in (A.4).

Analytic analysis of a steady-state Kalman filter estimator
The steady-state KF can be derived from the temporal continuous KF. If the system dynamic matrix F , the noise input mapping matrix G in (2), and the measurement matrix H are all constant matrices, and if the system noise and measurement noise are stationary ( Q and R are constant), the KF for a GNSS/INS/NHC integrated system may reach steady-state performance, leading to a constant covariance matrix P . In this condition, the Riccati equation for the continuous KF becomes an algebraic relation: In the steady-state condition, the rate at which the uncertainty increases (given by GQG T ) is balanced by the rate at which new information enters ( PH T R -1 HP ) and the dissipative effects of the system ( FP + PF T ). For a steadystate covariance matrix, the optimal filter is also time invariant, given by Taking the Laplace transform of this (neglecting initial conditions) yields so that x(s) is the optimal estimate of the state vector of the KF for the GNSS/INS/NHC integration. The term in brackets in the above equation is the transfer function representation of the steady-state KF. More details about the steady state can be found in the textbook [ (Maybeck , 1982), p. 273].

Navigation error propagation of the aided INS in the horizontal channel
The navigation error propagation analysis in the horizontal channel is similar to that for the vertical channel. According to assumption 1, the host vehicle is assumed to move at constant speed in a northward straight line. The north channel is the along-track direction, for which centimeter-level accuracy is sufficient and easy to fulfill when an accurate GNSS position is available. The east position component of the aided INS determines the track alignment measurement accuracy in the acrosstrack direction, which is critical for railway track geometry surveying (Chen et al. , 2015). Therefore, for the horizontal channel, we concentrate only on the east direction.
In addition to the position and velocity errors, the tilt error φ N is augmented into the error state vector, since a tilt error will introduce an acceleration component due to gravity with magnitude gφ N to be projected onto the horizontal axes (Titterton and Weston , 2004)). The attitude error about the yaw axis is also considered once the NHC is used. In this case the residual gyro bias about the x-axis δb gN and about the z-axis b gD , as well as the east component of the accelerometer bias δb aE , should also be augmented into the error state vector of the filter. where I 6 is the 6-by-6 identity matrix. δb gN is the north component of the gyro bias, modeled as a first-order Gauss-Markov process with correlation time T gb and driven by zero-mean white Gaussian noise w gbN ; δb gD is the vertical component of the gyro bias, modeled as a first-order Gauss-Markov process with correlation time T gb and driven by zero-mean white Gaussian noise w gbD . b aE is the east component of the accelerometer bias, modeled as a first-order Gauss-Markov process with correlation time T ab and driven by zero-mean white Gaussian noise w abE . w aE represents the accelerometer measurement noise along the east axis. w gN and w gD are the north and vertical components of the gyro measurement noise, respectively.
The measurement to be presented to the KF for the aided INS in the east direction should include the position difference between the INS and GNSS and the cross-track velocity difference between the INS and the NHC. Considering assumption 1, we have v D = 0 . Then, the velocity measurement Eq. (25) can be simplified as where n vE is the east component of n v in (12) and represents the measurement uncertainty of the across-track zero velocity. The GNSS position measurement equation can then be written as Therefore, the measurement equations for the aided INS in the east direction can be written as