Reliable Indoor Pseudolite Positioning Based on a Robust Estimation and Partial Ambiguity Resolution Method

The unscented Kalman filter (UKF) can effectively reduce the linearized model error and the dependence on initial coordinate values for indoor pseudolite (PL) positioning unlike the extended Kalman filter (EKF). However, PL observations are prone to various abnormalities because the indoor environment is usually complex. Standard UKF (SUKF) lacks resistance to frequent abnormal observations. This inadequacy brings difficulty in guaranteeing the accuracy and reliability of indoor PL positioning, especially for phase-based high-precision positioning. In this type of positioning, the ambiguity resolution (AR) will be difficult to achieve in the presence of abnormal observations. In this study, a robust UKF (RUKF) and partial AR (PAR) algorithm are introduced and applied in indoor PL positioning. First, the UKF is used for parameter estimation. Then, the anomaly recognition statistics and optimal ambiguity subset of PAR are constructed on the basis of the posterior residuals. The IGGIII scheme is adopted to weaken the influence of abnormal observation, and the PAR strategy is conducted in case of failure of the conventional PL-AR. The superiority of our proposed algorithm is validated using the measured indoor PL data for code-based differential PL (DPL) and phase-based real-time kinematic (RTK) positioning modes. Numerical results indicate that the positioning accuracy of RUKF-based indoor DPL is higher with a decimeter-level improvement compared that of the SUKF, especially in the presence of large gross errors. In terms of high-precision RTK positioning, RUKF can correctly identify centimeter-level anomalous observations and obtain a corresponding positioning accuracy improvement compared with the SUKF. When relatively large gross errors exist, the conventional method cannot easily realize PL-AR. By contrast, the combination of RUKF and the PAR algorithm can achieve PL-AR for the selected ambiguity subset successfully and can improve the positioning accuracy and reliability significantly. In summary, our proposed algorithm has certain resistance ability for abnormal observations. The indoor PL positioning of this algorithm outperforms that of the conventional method. Thus, the algorithm has some practical application value, especially for kinematic positioning.


Introduction
Pseudo-Satellite or Pseudolite, abbreviated as PL, is a transmitter deployed on the ground to transmit some kind of positioning signal, which usually transmits signals similar to navigation satellite system (GNSS) [1][2][3]. PL can be used as a supplement or enhancement for GNSS, especially in complex environments [4][5][6]. The PL system has been proven to have an independent positioning capability centimeter-level positioning precision in RTK under a good observation environment [27]. However, the indoor observation environment is often complex in practical application. The multipath effect of code observation and other interference factors may exist, and frequent observation abnormalities will negatively affect the accuracy and reliability of indoor PL positioning. Standard UKF (SUKF) lacks resistance to various abnormal situations [31,32]. Some scholars have put forward robust estimation theory to resist the influence of abnormal errors in measurement information [33,34]. The essence is to control the gross errors by constructing equivalent weights for weakening the influence of abnormal errors on the solution. This method reduces the contribution of anomalous observations to parameter estimation and can yield a relatively reliable result [35,36]. The recursive form of robust UKF (RUKF) is consistent with that of SUKF. The only difference between them is that the variance matrix of measurement noise in the observation model is replaced by the equivalent variance matrix [32,37]. Several familiar equivalent weight functions, such as IGG (I-III), Andrew, Tukey, and Huber methods, are usually adopted to calculate the equivalent weight matrix [38,39]. For indoor PL high-precision RTK, the RUKF will weaken the abnormal effects and improve the accuracy of float ambiguity solution when observation anomalies exist. However, all ambiguities cannot be guaranteed to be fixed successfully, especially under the influence of large gross errors. A partial AR (PAR) strategy is often adopted in the GNSS field to solve the aforementioned problem [40,41]. That is, if all ambiguity elements cannot be fixed reliably, then we can consider partial ambiguity elements that can be easily fixed. The PAR strategy is used to select the optimal ambiguity subset from all ambiguity elements and make a de-correlation for them. Then, the LAMBDA method is utilized to search the integer ambiguity solution for the selected ambiguity subset. Finally, the remaining ambiguity subset and coordinate solution are updated. In the GNSS field, the signal-to-noise ratio (SNR), elevation angle, posterior observation residual, the estimated variance matrix of float ambiguity, and the bootstrapping AR success rate or ratio value are usually used as the alternative information for selecting the optimal ambiguity subset [42]. Some scholars have successfully applied the robust Kalman filter to GNSS positioning and integrated navigation [43][44][45]. Most previous studies on PL positioning have focused on the outdoor simulated code-based meter-level positioning system and are thus not convincing [46][47][48]. Nearly no works have been performed on the RUKF-based indoor PL positioning, especially for the measured system and PAR-based high-precision RTK positioning.
In this study, we propose a reliable indoor PL positioning method based on the RUKF and PAR algorithm. First, the UKF is used for parameter estimation on the basis of posterior residuals. The standardized residues and anomaly recognition statistics are calculated. Then, the IGGIII scheme is used to construct the equivalent weight matrix for weakening the influence of abnormal observation. We further combine the RUKF with PAR strategy to obtain reliable PL positioning results with certain robustness by the anomaly recognition statistics information and select an ambiguity subset for PL-AR in the presence of large gross errors in the carrier phase observation. In Section 2, the specific method of our proposed method is provided. In Section 3, the superiority of the proposed method is verified using static and kinematic data for indoor DPL and RTK positioning. In Section 4, some conclusions are provided.

Indoor PL Positioning Model
For a single-frequency PL system, suppose that n + 1 PLs are observed synchronously on one epoch, 2n DD observation equations can be formed, namely, n DD code (p) and n DD carrier phase (φ) observation equations, which can be described as: ϕ ls qr = ρ ls qr + λN ls qr + ε ls qr p ls qr = ρ ls qr + e ls qr (1) where the subscripts q and r denote referenced and rover stations, respectively; the superscripts l and s denote referenced and non-referenced PLs, respectively; N ls qr denotes the PL DD ambiguity; the ε ls qr denotes DD carrier phase observation noise, and the e ls qr is DD code observation noise; For the indoor PL DD observations, the clock errors and phase center offset can be eliminated; ρ ls qr = (ρ s r − ρ l r ) − (ρ s q − ρ l q ) is the DD geometric distance between the receiver and PL. The mathematical model for the indoor PL relative positioning can be described as follows: where E[ • ] and D[ • ] denote the expectation and dispersion operators, respectively; a is the real-valued baseline vector; and b is the integer ambiguity vector. Γ denotes the design matrix that corresponds to a and Λ denotes the design matrix that corresponds to b. Q φφ denotes the covariance matrix that corresponds to the PL DD phase observation vectors, and Q pp denotes the covariance matrix that corresponds to the PL DD code observation vectors. The stochastic model for PL positioning generally uses an elevation-dependent weight model. The corresponding observation noise of the code and carrier phase can be pre-given approximate to GPS L1, and the noise can be obtained by prior modeling for PL observations. For the phase-based RTK model, if n + 1 PLs are synchronously observable, then n DD observation equation can be formed, and the contribution of code observation to high-precision RTK is generally ignored due to its poor accuracy. Thus, n + 3 parameters need to be estimated (n PL DD ambiguity and 3 coordinate parameters). Evidently, the number of single-epoch observations is less than the number of parameters to be estimated. For a PL system, multi-epoch static observation cannot be conducted to increase the number of observations due to strong correlation of static observations on successive epochs. Therefore, the coordinate should be regarded as an a priori value when PL ambiguity is solved in RTK. Thus, the KPI method is usually adopted for current PL systems.

Robust UKF
For a nonlinear system, the state and observation equation can be represented by: where k is the discrete time; x k is the state parameter to be estimated, in this study, it contains the coordinates and PL ambiguity parameters; y k is the measurement containing PL DD code and carrier phase observations; f (·) denotes the state function and is an identity matrix in this study; h(·) is a nonlinear observation function which can be seen as Equation (1); w k and v k are the corresponding Gaussian white noises which meet the following equation: where δ ij is the Kronecker − δ function; Q denotes the variance matrices of the process noise, R denotes the variance matrices of the measurement noise. The recursive form of RUKF is similar to that of SUKF. Their only difference is that the variance-covariance (VC) matrix of the measurement noise of the observation model is replaced by the equivalent variance matrix.
In this study, the calculation steps of RUKF algorithm are summarized as follows: Step 1: Initialization of the state parameters.
wherex 0 denotes the initial state parameters including the 3-D coordinate and PL DD ambiguity, and P 0 denotes the corresponding covariance matrices ofx 0 .
Step 2: Generation of sigma points. In this step, the sampling strategy based on Cholesky decomposition is usually utilized, for the n-dimensional variables x k−1 (with the mean and variance ofx k−1 and P k−1 ), the calculated vector χ k−1 that contains 2n + 1 sigma points is given as where D k−1 is the lower triangle matrix in the above Cholesky decomposition for P k−1 ; λ = α 2 (n + κ) − n is the scale factor, α and κ are constants, α determines the spread of the sigma points aroundx k−1 (usually set as 0 or 3 − n), and κ is a constant (usually set as 0). The corresponding weights of the produced sigma points are given as: where W n i denotes weight of mean value and W c i is weight of the covariance for the i-th sigma point; β is used to fuse prior information of random variables, and β= 2 is generally utilized for the Gauss distribution.
Step 3: Update the state.
wherex k|k−1 is the updated mean value of the state, and P k|k−1 denotes the corresponding covariance ofx k|k−1 .
Step 4: Update the measurement.
whereŷ k is the updated mean value of the predicted measurement, and P yy denotes the corresponding covariance ofŷ k .
Step 5: Update the SUKF filter where K k is the gain matrix;x k is the updated SUKF-based state parameter, and P k is the corresponding covariance matrix ofx k .
Step 6: Calculate the posterior observation residual In Step 5, we can obtain the estimated parameter dx (namely,x k above), then we can calculate the posterior observation residual v and corresponding VC matrix where B is the design matrix; in this study, B = Γ Λ Γ 0 as in Equation (2); l is a mis-closure vector of the DD PL observation, and C l is the VC matrix of l.
Step 7: Calculate the standardized residual The standardized residual for code and phase observations on the basis of v and C vv in Equation (11) can be described as where v k denotes the posterior observation residual of the k-th observation, C vv (k, k) is a scalar denoting the k-th diagonal element of C vv .
Step 8: Calculate the error discriminant statistics The constructed error discriminant statistics (∆v k ) on the basis of v k can be described as Step 9: Calculate the robust factor The calculated robust factor is used for adjusting the equivalent weight matrix. In this study, we adopt the IGGIII scheme. The three-segment function can be described as where k 0 and k 1 are constants, the reasonable ranges are generally 1.0 ≤ k 0 ≤ 2.0 and 2.5 ≤ k 0 ≤ 8.0.
Step 10: Update the RUKF filter We update the UKF filter solution by using the abovementioned robust factor.
where the meaning of P yy is the same as the expression in Equations (9) and (10). The subsequent filter updating is the same as Step 5.

PAR for PL Positioning
In the conventional indoor PL RTK processing, all PL DD ambiguities will be involved in the PL-AR using the LAMBDA method. The RUKF can improve the accuracy of the float solution to some extent, but the whole set of PL DD ambiguity cannot be fixed successfully. Thus, a PAR strategy for PL RTK positioning is introduced.
After conducting the RUKF, the parameters of the baseline vector and ambiguity float solution in Equation (2) denoteâ andb, and the corresponding VC matrix is composed of Qâ and Qb. The subset of the PL DD ambiguities can be described aŝ whereb is the entire set of estimated float ambiguity solutions;b 1 is the optimal subset, which can be easily fixed using the LAMBDA method; andb 2 is the remaining subset, which is difficult to fix due to certain errors in the corresponding DD phase observation. The selecting criterionb 1 is the core of the PAR for PL. In general, the SNR, elevation, posterior observation residual, the estimated variance information of float ambiguity, and the bootstrapping AR success rate or ratio value when conducting LAMBDA can be adopted to the identify the criterion of the optimal subset. In this study, we use the error discriminant statistics (∆v) in Equation (13) to select subsetb 1 .
We sort ∆v of all PLs in ascending order, which can be described as where ∆v k represents the elevation in k-th order. If the number of the subsetb 1 (m) is pre-given, then the program can easily obtain the subsetb 1 .
where N denotes the PL DD ambiguity float solution. If the subsetb 1 has been successfully fixed using LAMBDA, then the fixed solution denotes b 1 . Accordingly, we can update the remaining subsetb 2 , which can be described as If the subsetb 1 cannot be successfully fixed, then we continue to decrease the number of the subset b 1 (m) and repeat the above steps until m is less than 3. We then exit the program of PAR.

Data Processing
In this study, we propose a reliable indoor PL positioning method based on the RUKF and PAR. The relative positioning model and KPI method are adopted. If only code observation is used for positioning, we call it the DPL model, while if the high-precision carrier phase is used, we call it the RTK model. For the single-frequency PL system, the polynomial fitting method is used for cycle slip detection and reparation for carrier phase observations. The initial coordinate value is pre-given within a relatively accurate range when KPI is conducted. The popular LAMBDA method is used for PL-AR. Figure 1 shows the flowchart of our proposed method.

Observation Platform of the Indoor Positioning System
The PL positioning system is set up in a 10 m × 7 m × 4 m laboratory with five PLs mounted on the ceiling. The model of PL instrument is GSG-L1, which can generate an L1 carrier that is BPSK modulated with the C/A code and navigation signal. In this local coordinate, the origin point is located at the center of the laboratory. The coordinates of PLs were accurately determined in advance by a total station. As shown in Figure 2, the locations of PLs are indicated by red circles.

Observation Platform of the Indoor Positioning System
The PL positioning system is set up in a 10 m × 7 m × 4 m laboratory with five PLs mounted on the ceiling. The model of PL instrument is GSG-L1, which can generate an L1 carrier that is BPSK modulated with the C/A code and navigation signal. In this local coordinate, the origin point is located at the center of the laboratory. The coordinates of PLs were accurately determined in advance by a total station. As shown in Figure 2, the locations of PLs are indicated by red circles.

Observation Platform of the Indoor Positioning System
The PL positioning system is set up in a 10 m × 7 m × 4 m laboratory with five PLs mounted on the ceiling. The model of PL instrument is GSG-L1, which can generate an L1 carrier that is BPSK modulated with the C/A code and navigation signal. In this local coordinate, the origin point is located at the center of the laboratory. The coordinates of PLs were accurately determined in advance by a total station. As shown in Figure 2, the locations of PLs are indicated by red circles. The base and rover stations use Universal Software Radio Peripheral (USRP) as the frontend of the software receiver. The main function of USRP is to mix the L1 radio frequency signal to an intermediate frequency signal, and then complete the digitalization process. The following acquisition and tracking processes are all done by a self-developed software receiver, which brings great convenience to the research of indoor PL positioning systems. Figure 3 shows the experimental data acquisition scene. Each grid on the floor is a square of 0.6 m × 0.6 m. The base station is fixed on known points, and the rover station can move arbitrarily on a mobile car or along fixed rail. The base and rover stations use Universal Software Radio Peripheral (USRP) as the frontend of the software receiver. The main function of USRP is to mix the L1 radio frequency signal to an intermediate frequency signal, and then complete the digitalization process. The following acquisition and tracking processes are all done by a self-developed software receiver, which brings great convenience to the research of indoor PL positioning systems. Figure 3 shows the experimental data acquisition scene. Each grid on the floor is a square of 0.6 m × 0.6 m. The base station is fixed on known points, and the rover station can move arbitrarily on a mobile car or along fixed rail.

DPL Model
In this test, a static short baseline is used, and the observation conditions are relatively good. A total of approximately 11,000 observation epochs are used. Five PLs are synchronously observed during the entire observation period; thus, four DD code observations can be formed, we denote them as PL1-PL6, PL4-PL6, PL5-PL6 and PL8-PL6. The data processing method in Section 2.4 is used to conduct static robust DPL positioning for this test. The initial coordinates of UKF are provided accurately in the range of 0.1 m. A large gross error of 10 m is randomly added to some DD code observations on approximately 20 epochs to verify the resistance of the proposed algorithm to the gross errors of code observations and the small gross errors of the original observation data. Figure  4 shows the posterior residuals of four DD code observations and the corresponding sequence of error discriminant statistics ( i v Δ , denotes k as well) in Equation (13) during the entire observation period.

DPL Model
In this test, a static short baseline is used, and the observation conditions are relatively good. A total of approximately 11,000 observation epochs are used. Five PLs are synchronously observed during the entire observation period; thus, four DD code observations can be formed, we denote them as PL1-PL6, PL4-PL6, PL5-PL6 and PL8-PL6. The data processing method in Section 2.4 is used to conduct static robust DPL positioning for this test. The initial coordinates of UKF are provided accurately in the range of 0.1 m. A large gross error of 10 m is randomly added to some DD code observations on approximately 20 epochs to verify the resistance of the proposed algorithm to the gross errors of code observations and the small gross errors of the original observation data. Figure 4 shows the posterior residuals of four DD code observations and the corresponding sequence of error discriminant statistics (∆v i , denotes k as well) in Equation (13) during the entire observation period.
The observation posterior residuals can clearly reflect the anomalies of four PL DD code observations, especially large gross errors. A comparison of the error discriminant statistics (k) indicates that the corresponding values of k for some epochs with the 10 m gross error increase dramatically compared with those for other normal observation epochs. Apart from the epochs with artificial gross errors, many other epochs also show evident anomalies. This finding indicates that the error discrimination statistics can identify some small gross errors. In this experiment, k 0 and k 1 in Equation (14) are set to 2 and 8, respectively, when conducting the RUKF for indoor PL positioning. According to the absolute value of k, the observation environment can be divided into three situations as follows. In Situation #1, no abnormality is found in the observations, and the corresponding equivalent weight matrix does not need to be adjusted. In Situation #2, a certain small gross error exists in the observations, and the corresponding equivalent weight matrix needs weight reduction processing. In Situation #3, some large gross errors exist in the observatons, and the corresponding equivalent weight matrix needs weight abandonment processing. The specific method is presented in Equation (14). The observation posterior residuals can clearly reflect the anomalies of four PL DD code observations, especially large gross errors. A comparison of the error discriminant statistics ( k ) indicates that the corresponding values of k for some epochs with the 10 m gross error increase dramatically compared with those for other normal observation epochs. Apart from the epochs with artificial gross errors, many other epochs also show evident anomalies. This finding indicates that the error discrimination statistics can identify some small gross errors. In this experiment, 0 k and 1 k in Equation (14) are set to 2 and 8, respectively, when conducting the RUKF for indoor PL positioning. According to the absolute value of k , the observation environment can be divided into three situations as follows. In Situation #1, no abnormality is found in the observations, and the corresponding equivalent weight matrix does not need to be adjusted. In Situation #2, a certain small gross error exists in the observations, and the corresponding equivalent weight matrix needs weight reduction processing. In Situation #3, some large gross errors exist in the observatons, and the corresponding equivalent weight matrix needs weight abandonment processing. The specific method is presented in Equation (14). Table 1 shows the epoch number statistics of various DD PLs under the three situations mentioned above. The statistical results in Situation #1 indicate that most of the single DD code observations are normal without anomaly, but the epoch number of all the DD code observations that are simultaneously normal is relatively few at only 3694. The statistical results in Situation #2 show that most epochs, which account for approximately 63% of the total epochs, are affected by a minor gross error to varying extent. The statistical results in Situation #3 reveal that more than 300 epochs with a gross error of 10 m have relatively large gross errors in the original code observation.    Table 1 shows the epoch number statistics of various DD PLs under the three situations mentioned above. The statistical results in Situation #1 indicate that most of the single DD code observations are normal without anomaly, but the epoch number of all the DD code observations that are simultaneously normal is relatively few at only 3694. The statistical results in Situation #2 show that most epochs, which account for approximately 63% of the total epochs, are affected by a minor gross error to varying extent. The statistical results in Situation #3 reveal that more than 300 epochs with a gross error of 10 m have relatively large gross errors in the original code observation.  Figure 5 shows the positioning error sequence for SUKF and RUKF. The SUKF lacks resistance to observation anomalies, and its positioning results are worse than those of RUKF when gross errors exist for many epochs, especially large gross errors. This disadvantage of SUKF is evident. By contrast, the proposed RUKF algorithm can weaken the adverse effect of observation anomalies and improve the positioning accuracy. Figure 5 shows the positioning error sequence for SUKF and RUKF. The SUKF lacks resistance to observation anomalies, and its positioning results are worse than those of RUKF when gross errors exist for many epochs, especially large gross errors. This disadvantage of SUKF is evident. By contrast, the proposed RUKF algorithm can weaken the adverse effect of observation anomalies and improve the positioning accuracy.  Figure 6 shows the DPL average positioning error statistics (computed by the absolute value of every positioning errors) of SUKF/RUKF in Situations #2 and #3. A comparison of the statistical results shows that the RUKF performs better in Situations #2 and #3 than the SUKF, especially when large gross errors exist such as in Situation #3. In the case of small gross errors (Situation #2), RUKF has a decimeter-level improvement for the indoor PL positioning accuracy compared with SUKF because the RUKF uses a robust factor to adjust the equivalent weight matrix automatically. The RUKF also conducts weight reduction processing for anomalous observations, which can adaptively weaken its negative effect to some extent.  shows that the RUKF performs better in Situations #2 and #3 than the SUKF, especially when large gross errors exist such as in Situation #3. In the case of small gross errors (Situation #2), RUKF has a decimeter-level improvement for the indoor PL positioning accuracy compared with SUKF because the RUKF uses a robust factor to adjust the equivalent weight matrix automatically. The RUKF also conducts weight reduction processing for anomalous observations, which can adaptively weaken its negative effect to some extent.

RTK Model
The accuracy and stability of the DD carrier phase observation show much better performance than those of the DD code observation. The former is also less susceptible to interference and has smaller multipath effects under the indoor environment than the latter. However, if no robust strategies are applied under the occurrence of abnormal interference, then the final positioning accuracy and reliability will be affected; in particular, the PL-AR will be challenging. This

RTK Model
The accuracy and stability of the DD carrier phase observation show much better performance than those of the DD code observation. The former is also less susceptible to interference and has smaller multipath effects under the indoor environment than the latter. However, if no robust strategies are applied under the occurrence of abnormal interference, then the final positioning accuracy and reliability will be affected; in particular, the PL-AR will be challenging. This phenomenon is evident, especially for indoor kinematic positioning. In this study, the superiority of the proposed RUKF combined with the PAR algorithm is verified by static and kinematic tests in two cases. One case is that the PL ambiguity can be easily solved, and the other is that the the observation exists some abnormalities and PL ambiguity cannot be solved successfully using the conventional PL-AR method.

Static Test
In this experiment, one static short baseline is used, and the observation conditions are relatively good. A total of 270 observation epochs are used. Five PLs are observed synchronously during the entire observation period. Each epoch can form four carrier phase DD observations, denoted as PL4-PL1, PL5-PL1, PL6-PL1 and PL8-PL1. The PL4-PL1 observation for some epochs is added with a gross error of 5 cm without affecting the normal PL-AR to verify that the RUKF can identify and resist a certain small gross error for carrier phase observation. Figure 7 shows the posterior residuals and the corresponding error discrimination statistics (k) of the original PL4-PL1 observations, and the corresponding results after artificially adding a gross error of 5 cm. In Figure 7, the original carrier phase observations are relatively stable, most of the posterior residuals are less than 5 mm, and the corresponding k values are less than 2. Therefore, all epoch observations are not abnormal. After a small gross error of 5 cm is added, the corresponding epoch and posterior residual also exhibit a large jump. The corresponding k values are highly sensitive to detect the occurrence of the anomaly. When conducting KPI for the short baseline, the ambiguity is easily fixed as long as the given initial coordinates have high accuracy. The fix and hold mode is used to gain the ambiguity solution between successive epochs, and the coordinates for all epochs are fixed solutions in this test. Figure  8 shows the sequence of positioning errors for SUKF and RUKF. The SUKF-/RUKF-based positioning results have no differences in the epochs without gross error. However, for the epochs with a gross error of 5 cm (within the red ellipse in Figure 8), a centimeter-level jump occurs in the SUKF-based positioning results. By contrast, the RUKF weakens the influence of gross errors and does not show When conducting KPI for the short baseline, the ambiguity is easily fixed as long as the given initial coordinates have high accuracy. The fix and hold mode is used to gain the ambiguity solution between successive epochs, and the coordinates for all epochs are fixed solutions in this test. Figure 8 shows the sequence of positioning errors for SUKF and RUKF. The SUKF-/RUKF-based positioning results have no differences in the epochs without gross error. However, for the epochs with a gross error of 5 cm (within the red ellipse in Figure 8), a centimeter-level jump occurs in the SUKF-based positioning results. By contrast, the RUKF weakens the influence of gross errors and does not show large abnormality in its positioning results. The SUKS-based indoor PL-AR usually fails when gross errors exist in carrier phase observation. A gross error of 0.15 m is artificially added to the PL4-PL1 observations on the 100th epoch to verify the validity of our proposed RUKF and PAR strategy for indoor PL RTK positioning for the abovementioned static short baseline data. Three experimental schemes are adopted for data processing: Scheme #1 based on SUKF, Scheme #2 based on RUKF, and Scheme # 3 based on RUKF combined with the PAR algorithm. Figure 9 shows the positioning results for the three experimental schemes. For the 100th epoch, the SUKF-based PL-AR fails and has a jump of the positioning result of approximately 0.16 m due to the lack of resistance to abnormal observations. The RUKF cannot achieve the ambiguity fixed solutions as well, and shows a jump of the positioning result of approximately 0.12 m. Further comparison shows that SUKF needs 12 epochs to achieve ambiguity re-fixing, whereas RUKF needs five epochs. The reason is that the RUKF cannot realize PL-AR for the 100th epoch but still weakens the effect of adding gross errors to a certain extent. The accuracy of float ambiguity solution of RUKF is higher than that of SUKF and thus has a relatively faster filtering convergence speed for subsequent epochs. Scheme #3 adopts a PAR strategy based on RUKF. The PL-AR procedure does not consider all ambiguities for searching together, and an ambiguity subset that contains three other normal observations is selected for partial PL-AR. Therefore, the approach can achieve fixed solutions for all epochs and improves positioning accuracy. The SUKS-based indoor PL-AR usually fails when gross errors exist in carrier phase observation. A gross error of 0.15 m is artificially added to the PL4-PL1 observations on the 100th epoch to verify the validity of our proposed RUKF and PAR strategy for indoor PL RTK positioning for the abovementioned static short baseline data. Three experimental schemes are adopted for data processing: Scheme #1 based on SUKF, Scheme #2 based on RUKF, and Scheme # 3 based on RUKF combined with the PAR algorithm. Figure 9 shows the positioning results for the three experimental schemes. For the 100th epoch, the SUKF-based PL-AR fails and has a jump of the positioning result of approximately 0.16 m due to the lack of resistance to abnormal observations. The RUKF cannot achieve the ambiguity fixed solutions as well, and shows a jump of the positioning result of approximately 0.12 m. Further comparison shows that SUKF needs 12 epochs to achieve ambiguity re-fixing, whereas RUKF needs five epochs. The reason is that the RUKF cannot realize PL-AR for the 100th epoch but still weakens the effect of adding gross errors to a certain extent. The accuracy of float ambiguity solution of RUKF is higher than that of SUKF and thus has a relatively faster filtering convergence speed for subsequent epochs. Scheme #3 adopts a PAR strategy based on RUKF. The PL-AR procedure does not consider all ambiguities for searching together, and an ambiguity subset that contains three other normal observations is selected for partial PL-AR. Therefore, the approach can achieve fixed solutions for all epochs and improves positioning accuracy.

Kinematic Test
In this test, a fixed rail is used, and a mobile car with a PL receiver antenna travels from one end to the other. A total of 270 observation epochs are used. Five PLs are observed synchronously during the entire observation period. Each epoch can form four carrier phase DD observations. The advantage of our proposed method is verified in the case of successful and failed PL-AR with conventional methods.
A gross error of 0.1 m is added to one carrier phase observation for some epochs, and the SUKF and RUKF are used for kinematic positioning. The initial coordinates given for the first epoch are relatively accurate, the fix and hold mode is adopted to gain the ambiguity solution between successive epochs, and the coordinates of all epochs are fixed solutions in this experiment. Figure 10 shows the positioning results of two different schemes. The SUKF-based X-directional positioning results show several centimeter-level jumps for the epochs with gross errors (within the red ellipse in Figure 10), whereas RUKF can effectively weaken the effect of adding gross errors and has no large fluctuation in its positioning results.

Kinematic Test
In this test, a fixed rail is used, and a mobile car with a PL receiver antenna travels from one end to the other. A total of 270 observation epochs are used. Five PLs are observed synchronously during the entire observation period. Each epoch can form four carrier phase DD observations. The advantage of our proposed method is verified in the case of successful and failed PL-AR with conventional methods.
A gross error of 0.1 m is added to one carrier phase observation for some epochs, and the SUKF and RUKF are used for kinematic positioning. The initial coordinates given for the first epoch are relatively accurate, the fix and hold mode is adopted to gain the ambiguity solution between successive epochs, and the coordinates of all epochs are fixed solutions in this experiment. Figure 10 shows the positioning results of two different schemes. The SUKF-based X-directional positioning results show several centimeter-level jumps for the epochs with gross errors (within the red ellipse in Figure 10), whereas RUKF can effectively weaken the effect of adding gross errors and has no large fluctuation in its positioning results.
In kinematic positioning of a moving rover station, the initial coordinate value of the current epoch is usually based on the coordinates of the previous epoch. A large deviation in the positioning result of one epoch negatively affects the subsequent initialization of coordinates and thus leads to filtering divergence and positioning failure. Therefore, the proposed RUKF and PAR algorithm is more necessary than the static positioning. In this experiment, a gross error of 0.15 m is added to one carrier phase observation for the 50th epoch. Two different data processing schemes, namely, SUKF and RUKF and PAR, are adopted for kinematic positioning. The single-epoch PL-AR model is used during the successive epochs. Figure 11 shows the planar trajectory derived from the kinematic positioning results of the two different schemes. For the 50th epoch, the SUKF-based PL-AR fails and has a deviation of several decimeters for the positioning result, which is transmitted to the subsequent initialization of coordinates. The subsequent filtering of epochs cannot converge and continues to diverge along with the corresponding positioning results due to the large deviation. Figure 11 shows that the plane trajectory correspondingly deviates from the true straight rail. When the RUKF and PAR algorithm is adopted, the effect of adding the gross error is weakened to a certain extent, and the ambiguity subset containing three other normal observations is selected and successfully realizes the partial PL-AR. The positioning accuracy based on our proposed method is better by 2-3 cm than the planar positioning trajectory and the true straight rail. In kinematic positioning of a moving rover station, the initial coordinate value of the current epoch is usually based on the coordinates of the previous epoch. A large deviation in the positioning result of one epoch negatively affects the subsequent initialization of coordinates and thus leads to filtering divergence and positioning failure. Therefore, the proposed RUKF and PAR algorithm is more necessary than the static positioning. In this experiment, a gross error of 0.15 m is added to one carrier phase observation for the 50th epoch. Two different data processing schemes, namely, SUKF and RUKF and PAR, are adopted for kinematic positioning. The single-epoch PL-AR model is used during the successive epochs. Figure 11 shows the planar trajectory derived from the kinematic positioning results of the two different schemes. For the 50th epoch, the SUKF-based PL-AR fails and has a deviation of several decimeters for the positioning result, which is transmitted to the subsequent initialization of coordinates. The subsequent filtering of epochs cannot converge and continues to diverge along with the corresponding positioning results due to the large deviation. Figure 11 shows that the plane trajectory correspondingly deviates from the true straight rail. When the RUKF and PAR algorithm is adopted, the effect of adding the gross error is weakened to a certain extent, and the ambiguity subset containing three other normal observations is selected and successfully realizes the partial PL-AR. The positioning accuracy based on our proposed method is better by 2-3 cm than the planar positioning trajectory and the true straight rail.

Conclusions
Our observations in this study are prone to abnormalities due to the complex environment of indoor PL positioning. The traditional parameter estimation method lacks the effective ability of anomaly recognition and anti-interference, and the conventional PL-AR becomes challenging. This study presents a reliable indoor PL positioning method based on the RUKF and PAR strategy. The error discriminant statistic is constructed using posterior residuals of the observations, and the equivalent weight matrix of observation is adaptively adjusted by the calculated robust factor. The

Conclusions
Our observations in this study are prone to abnormalities due to the complex environment of indoor PL positioning. The traditional parameter estimation method lacks the effective ability of anomaly recognition and anti-interference, and the conventional PL-AR becomes challenging. This study presents a reliable indoor PL positioning method based on the RUKF and PAR strategy. The error discriminant statistic is constructed using posterior residuals of the observations, and the equivalent weight matrix of observation is adaptively adjusted by the calculated robust factor. The optimal ambiguity subset is selected to conduct the partial PL-AR on the basis of the error discriminant statistics. The superiority of the proposed algorithm for indoor code-based DPL and phase-based RTK positioning is verified using static and kinematic observation data. Some valuable conclusions are obtained as follows.

1.
Compared with the SUKF algorithm, the RUKF algorithm can effectively weaken the anomalous effect of PL code observations and improve the accuracy and reliability of indoor DPL positioning, especially when certain large gross errors exist.

2.
RUKF can identify small gross errors (centimeter-level) of PL carrier observations and achieve the corresponding indoor RTK positioning accuracy of fixed solutions at a centimeter-level improvement.

3.
Compared with SUKF, RUKF can improve the accuracy of ambiguity float solution and the re-convergence speed when the carrier phase observation has relatively large gross errors. However, RUKF cannot achieve PL-AR successfully. The proposed RUKF combined with PAR strategy can achieve partial PL-AR for the selected ambiguity subset and obtain an accurate fixed solution.
The advantages of our proposed algorithm are important for indoor PL kinematic positioning.
Currently, most of the PL systems are still only capable of high-positioning adoption of the KPI, which brings some inconveniences in application. In our future work, we will form a combination for indoor PL and ultra-wideband (UWB) positioning, and it is expect that the indoor UWB-assisted PL positioning can avoid the disadvantage of KPI and have a better positioning performance.
Author Contributions: X.L., G.H., and P.Z. conceived the study and designed the review. X.L. and P.Z. conducted the literature review. X.L. wrote the first draft. P.Z. performed the data collection. Q.Z. provided critical comments and fund support for the study. All the authors contributed to the writing and reviewing of the manuscript.