- Email: [email protected]
Switching position and range-domain carrier-smoothing-code filtering for GNSS positioning in harsh environments with intermittent satellite deficiencies
Available online at www.sciencedirect.com
Journal of the Franklin Institute 356 (2019) 4928–4947 www.elsevier.com/locate/jfranklin
Switching position and range-domain carrier-smoothing-code filtering for GNSS positioning in harsh environments with intermittent satellite deficiencies Guobin Chang a,c,∗, Tianhe Xu b,c, Chao Chen a, Bing Ji d, Shengquan Li e a School
of Environmental Science and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China b Institute of Space Science, Shandong University, Weihai 264209, China c State Key Laboratory of Geo-information Engineering, Xi’an Research Institute of Surveying and Mapping, Xi’an 710054, China d Department of Navigation Engineering, Naval University of Engineering, Wuhan 430033, China e Acoustic Science and Technology Laboratory, College of Underwater Acoustic Engineering, Harbin Engineering University, Harbin 150001, China Received 10 May 2018; received in revised form 19 December 2018; accepted 2 April 2019 Available online 11 April 2019
Abstract Carrier-smoothing-code filtering (CSCF) is widely used in GNSS signal processing to combine code pseudoranges and carrier phases. Position-domain (PD) CSCF is generally more accurate and less sensitive to visible satellite changes than range-domain (RD) CSCF. However, PD-CSCF necessitates at least four visible satellites. Intermittent satellite deficiency with less than four visible satellites is not uncommon in harsh environments like urban canyons. At such deficiency epochs, the PD-CSCF convergence has to break off. This study aims to bridge intermittent deficiency epochs in PD-CSCF without introducing any external information. The proposed solution is called switching RD and PD CSCF in which PD filter is replaced by RD filter at deficiency epochs. Besides detailing the seamless switching algorithms from PD/RD to RD/PD filters, a global RD filter, different from the conventional one, is developed to preserve the correlations among smoothed pseudoranges corresponding to different satellites. Compared to the conventional PD and RD CSCF algorithms, smoother results can be expected from the proposed switching filter, especially after the intermittent deficiency epochs. Experiments are ∗ Corresponding author at: School of Environmental Science and Spatial Informatics, China University of Mining and Technology. E-mail address: [email protected] (G. Chang).
https://doi.org/10.1016/j.jfranklin.2019.04.005 0016-0032/© 2019 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4929
conducted using real BDS signals. Cases with different kinds of deficiency are considered. Superiority of the proposed method is clearly observed from the results. © 2019 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.
1. Introduction Global Navigation Satellite Systems (GNSS) technology is widely used in broad engineering applications relating navigation, guidance and control. Basically, two types of signals are normally used in GNSS, namely code pseudoranges and carrier phases, though Doppler signals are also incorporated in some applications [1]. Compared with code pseudorange, carrier phase is with much higher accuracy but ambiguous with unknown integer cycles. In order to make full use of the high-accuracy carrier phase signals, two approaches can be used to address the ambiguity problem, namely performing integer ambiguity resolution (AR) including both estimation and validation and smoothing the code pseudoranges with carrier phases [2]. The latter, called the carrier-smoothing-code filtering (CSCF) technique, is our focus in this work. If we include integer cycles as parameters to be estimated as float values without fixing them as integers, the corresponding approach can roughly be regarded as a CSCF technique, though it is called carrier adjusted code (CAC) in [3]. In the CAC approach when applied in dynamical positioning, the integer cycles function as kind of media to carry historical information forward to current and future epochs. Of course, AR and CSCF approaches can be merged, e.g., before integer cycles can be fixed with sufficiently high success rate, CSCF can be used in the transient stage, or in cases of multiple-frequency signals, one can perform partial resolution of the integer ambiguities, namely linear combinations of the original integer ambiguities, with the remaining part dealt with using CSCF technique [4]. In CSCF, what one does is nothing but a combination of the code pseudoranges and the time differenced or epoch differenced carrier phases. The aim of the time difference operation is to cancel integer ambiguities as long as cycle slips are absent or repaired. Cycle slips detection and correction is a prerequisite for almost all applications involving carrier phases, not only for the CSCF studied in this work, but also for AR. Cycle slip detection/correction can and also have to be performed at different stages of the whole process of data processing, it is far from a trivial task, especially when performed in real time [5]. However, to be more focused, the cycle slip repairing issue is not to be explored in this study and we simple assume the absence of any cycle slips in the used carrier phase data. Time differenced carrier phases, carrying only relative or incremental information of the position, can be solely relied on in determining displacements, velocities or accelerations. However, they have to be combined with code pseudoranges to obtain non-divergent absolute positioning solutions, similarly as dead-reckon positioning techniques like inertial navigation system. The combination of code pseudoranges and time differenced carrier phases in the CSCF can be carried out generally in two different ways, namely in range domain (RD) and in position domain (PD). In RD-CSCF, signals are processed separately for different satellites [2,6–8]. However, in PD-CSCF, all available measurements, namely the code pseudoranges and the time differenced carrier phases from all visible satellites (considering the minimum/mask elevation angle), are processed altogether to estimate the position coordinates and of course the receiver clock error. This is nothing but exactly a process of information fusion. Although other fusion strategies or some extra tunings
4930
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
may be employed [9,10], the fusion is often optimally performed in the least-squares sense [11–13]. When used in static positioning, e.g., at the reference stations of ground or satellite based augmentation systems, this filtering degenerates to a recursive least-squares estimation. It has been shown, both theoretically and experimentally, PD-CSCF can be superior to RDCSCF in general in terms of accuracy and insensitivity to visible satellite changes [10,13]. This is not any surprise when considered in the framework of information fusion or statistical signal processing. In PD-CSCF, besides the information propagation along time line, further information sharing across different satellites/channels is also permitted. Satellite redundancy, not merely sufficiency, namely with more than 4 visible satellites, can be exploited fully by PD-CSCF which RD-CSCF fails to do. The more visible satellites there are, the more obvious of the advantage of the PD-CSCF will be. Note that the following needs to be emphasized to avoid misunderstanding: with exactly 4 satellites visible, namely without redundancy, PDCSCF and RD-CSCF would not be equivalent to each other either, in general, or to be more specific, the former would still outperform the latter. The more accurate PD-CSCF is also more restrictive. At least four satellites have to be visible in order to perform the PD-CSCF. As the single-channel RD-CSCF is carried out separately for each visible satellite, there is no such restriction on the number of visible satellites, though positioning using the output of the RD-CSCF, namely the smoothed pseudoranges, can be possible only with at least four visible satellites. For scenarios with intermittent satellite deficiencies, RD-CSCF would outperform PD-CSCF in a sense as explained in the following. Although neither can provide positioning service at the deficiency epochs, RD-CSCF can still continuously pass the historical information buried in the still visible satellites on to the smoothed pseudoranges at these deficiency epochs and also on to those at future sufficiency epochs; however, PD-CSCF can NOT ensure such continuous information passing. Satellite deficiency at any epoch would break off the converging process of the PD-CSCF. The historical information before the deficiency epoch has to be entirely abandoned, and the convergence has to start over again when satellite sufficiency restores. Frequent deficiencies would break the filtering process into pieces in which smoothing effects could not be obvious owing to the short converging periods. This of course would result in poor positioning performances. Scenarios of satellite deficiency are not uncommon in some harsh environments with blocking/jamming or multipath, such as urban canyons. In the past, this deficiency is mitigated mainly by fusing additional information [14,15]. The additional information can be those provided by external sensors and/or infrastructures. Exploiting the model of the environment, the vehicle or the equipment/sensor, such as the road maps, the vehicle kinematics models, and the receiver clock models can also provide information to permit a solution during satellite deficiencies. Of course these efforts are meaningful, however at a price, e.g., cost and/or complexity would be increased. There are some situations in which short period of positioning unavailability can be tolerated, e.g., when a car drives under a viaduct or passes by a tall building. In these situations, cumbersome job of adding external sensors or exploiting models may not be cost effective. In this work it is tried to generalize PD-CSCF to the case with intermittent satellite deficiencies and without any external information. The solution, seemingly straightforward, can be using a switching PD and RD CSCF approach, namely using PD-CSCF for satellitesufficiency epochs and RD-CSCF for satellite-deficiency epochs. Then it is worthwhile to conduct a seamless transition from PD-CSCF to RD-CSCF and of course also a seamless opposite transition. By saying seamless, we mean that with the transition from PD to RD, all relevant statistical information in the coordinates/clock-error estimate in the PD-CSCF would
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4931
be propagated into the smoothed pseudoranges in the RD-CSCF and vice versa. To our best knowledge, this switching strategy has not been studied in the literature. Our aim to conduct this study is to present the detailed model and algorithm for the switching PD and RD CSCF for GNSS-only applications in cases with intermittent satellite deficiencies. While the PD algorithm for epochs of satellite sufficiency in the proposed switching PD and RD CSCF method is exactly the same to the conventional PD-CSCF, the RD algorithm for epochs of satellite deficiency deviates from the traditional RD-CSCF algorithm. At the transition epoch from satellite sufficiency to deficiency, the same set of statistical information at previous epoch flows into any one of the smoothed pseudoranges at the current epoch. Apparently, the thus obtained smoothed pseudoranges are no longer independent in contrary to the conventional channel-wise RD-CSCF. At the subsequent epochs with satellite deficiency, if we conduct the conventional single-channel RD filtering for each visible satellite separately as in the conventional RD-CSCF, the solution would not be statistically optimal owing to ignoring the correlations among different channels. So instead, in the proposed switching algorithm, the RD processing at deficiency epochs for all visible satellites are performed as a whole rather than separately. This is statistically correct and brings in global optimality. At deficiency epochs, both the switching method and the conventional RD-only method are conducting RD filtering, however, each one of the smoothed pseudoranges provided by the latter only contain the historical information of the corresponding satellite while the former also contain the historical information of other, previously visible satellites. This would produce smoother results enabling higher accuracy and faster convergence. This is exactly what we are pursuing by conducting this study. The single-frequency point positioning is employed as an example in both the method presentation and the experiment study. Assume all error terms except clock error and noises have been corrected using models. This simplifies the model with most error terms absent and makes us more focused on the key element, namely the switching PD and RD CSCF algorithm development and validation. This is of course without loss of generality; the proposed method can be readily extended to other cases, e.g., relative (differential) positioning and/or multiple frequencies. For relative positioning with short to mediate baselines, the problem becomes even easier with most common error terms being cancelled or reduced largely [8]; and for multiple frequencies, linear combinations of raw measurements may be used as the inputs of the CSCF techniques including the one proposed in this work [4,6,7]. Of course, the combination and the CSCF can be performed in one step where global optimality can be achieved by selecting the combining and filtering coefficients simultaneously. The rest of the paper is organized as follows. In Section 2, the conventional PD and RD CSCF algorithms are presented in order to introduce terminology and to make the proposed method in the context. The main algorithm is presented in Section 3. Experiments using real BDS data are presented in Section 4 and the results are analyzed to check the performance of the proposed method. Concluding remarks are given in section V also with some potential topics of future study suggested.
2. Range-domain and position domain carrier-smoothing-code filters Kinematic point positioning is employed here. Error terms are assumed to be absent or removed/corrected using empirical models, such as satellite orbit and clock errors, tropospheric and ionospheric delay errors.
4932
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
2.1. Measurement model The code pseudorange and carrier phase (in meters) measurement equation would then be in the following simple forms [3]: ρkj = rkj + tk + ekj
(1)
ϕkj = rkj + tk + λnk + εkj
(2)
In the above, subscript k and superscript j denote the indices of epoch and the satellite, respectively; ρ and ϕ denote the code pseudoranges and carrier phases, respectively and e and ε denote the corresponding noises including all un-modeled errors, such as hardware noises and multipath errors, satellite orbit and clock errors and residual tropospheric/ionospheric delay errors; λ and n denote the carrier frequency and the unknown integer cycles, respectively; t denotes the receiver clock error multiplied by the light speed; r denotes the geometric range whose expression is as follows: j 2 2 2 xk − xk + ykj − yk + zkj − zk rkj = (3) where x/y/z denotes the Cartesian coordinate triad; the x/y/z with superscript j denote the known coordinates of the satellite and those without any superscript denote the coordinates of the rover receiver which, together with tk , are the parameters to be estimated. In most CSCF approaches, time differencing carrier phases is firstly conducted to remove the unknown but constant integer cycles; so the following measurement equation is used instead of Eq. (2): j δϕkj = δrkj + δtk + εkj − εk−1
(4)
Of course Eq. (4) relies on the absence of cycle slips. Preprocessing of detecting and repairing the potential cycle slips is needed, which is out of the range of this study. 2.2. Range-domain carrier-smoothing-code filters In RD-CSCF, the range plus receiver clock error, namely, pkj = rkj + tk
(5)
called (true) pseudorange, is the parameter firstly to be estimated, though the coordinates of the rover station are of final interest. The first and most widely used version of this kind of approach is the Hatch filter, which is a first-order infinite impulse response filter with fixed or varying but predetermined filtering coefficients [2]. However, in this work, a stochastic filtering methodology based on the least-squares criterion, rather than the digital filtering in the original Hatch filter, is followed. There are no essential differences between the stochastic and the Hatch filters; we are only replacing the predetermined filtering coefficient in the latter with a variant one adaptively tuned in the former. This variant coefficient is dynamically calculated according to the variances of the involved variables. The filtering algorithm is presented for the initialization and the recursion stages separately for convenience. For the initialization stage, the output, namely the smoothed pseudorange, as kind of estimate of the (true) pseudorange, denoted as pˆ kj , is assigned directly as the code pseudorange. Of course the variances of the smoothed pseudoranges are assigned as those of the code pseudoranges. Then for the case at the initialization epoch, the covariance defined in the following is zero: ckj = cov pˆ kj , ϕkj = 0
(6)
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4933
This is clearly reasonable, because there is not any carrier phase signals involved in calculating pˆ kj at the initialization epoch. Note that the epoch k does not necessarily equate to zero or one. Any epoch at which a satellite starts to be visible is an initialization epoch for this satellite. Any epoch at which cycle slip is detected but failed to be corrected (considering the success rate) is also viewed as an initialization epoch. The algorithm for the recursion is given in the following. For each epoch, there are two kinds of measurement for estimating the pseudorange. The first is that propagated from the previous epoch using the time differenced carrier phases, namely: j j j p¯ kj = pˆ k−1 +δϕkj = pkj + ξk−1 + εkj − εk−1 = pkj + ςkj j ξk−1
(7)
j pˆ k−1
2 where denotes the error in with the corresponding variance denoted as σˆ k−1 , j . The variance of the error in the above measurement is as follows: j j j j + var εkj + var εk−1 − 2cov ξk−1 σ¯ k,2 j = var ςkj = var ξk−1 , εk−1 j 2 2 2 = σˆ k−1 , j + σϕ,k, j + σϕ,k−1, j − 2ck−1
(8)
j 2 2 2 where σˆ k−1 , j and ck−1 are calculated at the previous epoch and σϕ,k, j /σϕ,k−1, j are the variances of the carrier phase measurements which are assigned empirically and based on elevation angle. The second kind of measurement is the raw code pseudorange as shown in Eq. (1). Combine the measurements in Eqs. (1) and (7) optimally in the least-squares sense as follows:
j −2 pˆ kj = σˆ k,2 j σ¯ k,−2j p¯ kj +σρ,k, j ρk
(9)
with
−2 −1 σˆ k,2 j = var pˆ kj = σ¯ k,−2j +σρ,k, j
(10)
2 ck = cov pˆ kj , ϕkj = cov σˆ k,2 j σ¯ k,−2j ϕkj , ϕkj = σˆ k,2 j σ¯ k,−2j σϕ,k, j
(11)
2 where σρ,k, j is the variance of the code pseudorange measurements which is also assigned empirically and based on elevation angle. Note that besides the variances in Eq. (10), the cross covariance calculated in Eq. (11) is necessary for the next epoch as is apparent in Eq. (8). At any epoch, with the smoothed pseudoranges available from at least four satellites, the coordinates of the rover station can be determined together with the receiver clock error, in the same way as the standard code pseudorange positioning. Note that all error terms common to code pseudoranges and carrier phases, e.g., the satellite orbit and clock errors and the residual tropospheric delay errors, do not matter in the determination of the smoothed pseudoranges, because all these errors are absorbed in the pseudoranges. It is the coordinate determination in which these errors enter as factors; however, they are assumed small after the model correction in this work.
2.3. Position-domain carrier-smoothing-code filters Contrary to the two-step strategy in the RD-CSCF algorithm presented in the above, namely first smoothing pseudoranges and then determining coordinates/clock-error, the PD-CSCF employs a one-step strategy, namely, determining the coordinates/clock-error directly using all available code pseudoranges and time differenced carrier phases. Let βk = [xk yk zk tk ]T be the parameter vector. Similarly to the RD-CSCF, at the initialization epoch, the parameter vector
4934
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
is estimated solely using code pseudoranges, and the cross covariance matrix between the parameter estimates and the carrier phase measurement from the jth satellite at the same initialization epoch, denoted as ukj = cov[βˆk ,ϕkj ], is a 4 × 1 zero vector. Note that visible satellite set changes will not trigger re-initialization as long as there are still at least four satellites visible. This is exactly the characteristic that PD-CSCF is not sensitive to satellite changes as stated above. However, at any epoch with satellite sufficiency restoring from previous deficiency, the PD-CSCF will re-initialize and hence this epoch is also an initialization epoch for the conventional PD-CSCF. For detailed description of PD-CSCF see e.g., [10,13]. Substitute Eq. (3) into Eq. (1); perform linearization at some approximate values of the coordinate; and we have the following measurement equation: lkj = aTk, j βk + ekj
(12)
where l and a denote the observation and design coefficient, respectively, with the superscript and subscript defined as before. Note that in Eq. (12), we directly treat the parameter vector βk as the argument rather than the correction or delta vector, for convenient presentation of the method. Either the coordinates from the standard positioning using only code pseudoranges or from the solution at the previous epoch can be used as the approximate values for the linearization, namely calculating the a vectors in Eq. (12). Stack equations as in Eq. (12) for all visible satellites, we have the vector measurement equation: lk = Ak βk + ek
(13)
Then the least-squares estimation is as follows: −1 T −1 β˜k = ATk Q−1 Ak Qee,k lk ee,k Ak
(14)
T −1 −1 Qβ˜ β,k ˜ = Ak Qee,k Ak
(15)
This represents the first piece of information for the parameter. The second piece is derived as follows. Linearize Eq. (4) with respect to time differenced vector δβk as follows: j j T ykj = ak, j δβk + εk − εk−1
(16)
where y denote the observation; and a is the same as defined in Eq. (12). Stack all available equations as in Eq. (16), we have the vector measurement equation as follows: yk = Bk δβk + εk − εk−1
(17)
Note here though the coefficient a in Eq. (16) is the same as that in Eq. (12), the coefficient matrix B in Eq. (17) does not necessarily be the same as that in Eq. (13). The visible satellite set represented in Eq. (16) may not be the same as that in Eq. (12), because a satellite represented in the former has to be continuously tracked at both the previous and the current epochs, while that in the latter only needs to be tracked at the current epoch. The least-squares estimation for Eq. (17) is as follows: −1 δ β¯k = BTk Qεε,k + Qεε,k−1 Bk −1 BTk Qεε,k + Qεε,k−1 −1 yk = Fk yk (18) T −1 −1 Qδβδ Bk ¯ β,k ¯ = Bk Qε ε ,k + Qε ε ,k−1
(19)
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
Let fk,j denote the jth column of Fk . Then we have the following cross covariance: ⎡ ⎤ j j ⎦ Gk = cov βˆk−1 , δ β¯k = cov⎣βˆk−1 , − fk, j εk−1 =− uk−1 fk,T j j∈ 0
4935
(20)
j∈ 0
In the above 0 denotes the index set of all the visible satellites from which the time differenced carrier phase signals are both available at the (k − 1)th and the kth epochs. Or in other words, a satellite whose index is in 0 has to be continuously tracked at the (k − 2)th, the (k − 1)th and the kth epochs. Then we have the propagated estimates for the parameter at the kth epoch as follows: β¯k = δ β¯k + βˆk−1
(21)
Qβ¯ β,k + GkT + Gk ¯ = Qδ βδ ¯ β,k ¯ + Qβˆ β,k−1 ˆ
(22)
This represents the second piece of information for estimating the parameter at the current epoch. It is clear that these two pieces of information are independent. So by simply fusing them optimally in the least-squares sense, we have the following expression for the final parameter estimates at the current epoch:
−1 −1 −1 Q Qβˆ β,k = + Q (23) ˆ ¯ ˜ β¯ β,k β˜ β,k
−1 ¯ Qβ−1 βˆk = Qβˆ β,k β k ˆ ¯ β,k ¯ βk + Qβ˜ β,k ˜
(24)
j j −1 −1 2 ukj = cov βˆk , ϕkj = cov Qβˆ β,k Q f ϕ , ϕ k, j ˆ ˆ Qβ¯ β,k ¯ ¯ f k, j σϕ,k, j k k = Qβˆ β,k β¯ β,k
(25)
Again, similar to the RD-CSCF, besides the covariance matrix in (23), the cross covariance (vector) in (25) is also necessary which will be used in the estimation at the next epoch. This cross covariances exactly represents the specific consideration of the time correlations in the time difference carrier phase signals. The PD-CSCF algorithm presented here is exactly the same to that derived in [12,13], though with a slightly different presentation. The following two are noted concerning the PD-CSCF. First, the second piece of information, namely the estimate in Eq. (21) may be called a prediction in some references. However, this cannot be a true prediction as that in the conventional Kalman filter. As is clear from the above derivation, this estimate, namely the so called prediction, can only be available at the current epoch, the same to the estimate in Eq. (14), because the carrier phase signals at the current epoch are involved in constructing the time differenced carrier phase signals. This is contrary to the conventional Kalman filter in which the prediction or time update is carried out at the (k − 1)th epoch [16,17]. So in this sense, there is no prediction of any physical meaning at all in the PD-CSCF. Second, from the presentation given here, both the number of the code pseudoranges and that of the time differenced carrier phases have to be at least four. However, this is an excessive demand. In fact either of the two numbers being larger than or equal to four can be sufficient for the existence of a parameter estimate at the current epoch. Or even less stringently, the PD-CSCF algorithm will be applicable as long as measurements from at least four different satellites are available; no matter the measurements are the code pseudoranges or the time differenced carrier phases. This is clearer when the PD-CSCF algorithm is presented in the phase-connected-code form [11]. Finally after the presentation, we want to stress once
4936
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
more that the PD-CSCF, no matter in what form, necessitates that at least four satellites have to be visible above the predetermined minimum/mask elevation angle. The CSCF algorithm for cases with intermittent satellite deficiencies, which is the main contribution of this work, is developed in detail in the next section.
3. Switching range and position domain carrier-smoothing-code filter In this section, the switching PD and RD CSCF algorithm is presented. As stated before, the RD algorithm in the proposed switching method at deficiency epochs deviates from the conventional single-channel RD algorithm presented in the previous section; so the RD algorithm here, essentially a global one, is also detailed in this section. Without loss of generality, the switching algorithm is initialized for a satellite-sufficiency epoch as the conventional PDCSCF presented in the last part of the previous section. As a result, in the following, the switching algorithm is presented in three parts, namely for the transition from the PD to the global RD filters, for the global RD filtering at the intermittent satellite-deficiency epochs, and for the transition from the global RD to the PD filters, respectively.
3.1. Switching from position-domain to global range-domain filters Once the number of visible satellites becomes deficient, namely with sufficient visible satellites at the immediately previous epoch, the transition from the PD to the global RD CSCF is triggered. The input to this transition is exactly the same to normal PD algorithm, namely the parameter estimate, the associate covariance matrix and the cross covariance vectors between the parameter estimate and the involved carrier phases. At this transition epoch, the parameters to be estimated changed from the coordinates/clock-error to the pseudoranges. Stack all the pseudoranges corresponding to the still visible satellites as the pseudorange vector, denoted as pk . Note that for 2 or 3 visible satellites, this vector are indeed a vector, with 2 and 3 elements, respectively. However, for the case of only one visible satellite, this vector degenerates to a scalar. Similarly stack code pseudorange and carrier phase measurements for all satellite still visible at the current epoch, in the same order as in pk , and denote them as the code pseudorange vector ρk and the carrier phase vector ϕk , respectively. For estimating the pseudorange vector, the code pseudoranges provide a piece of information. The time differenced carrier phases provide a piece of information for the relative or incremental pseudoranges in the interval from the previous to the current epoch. With the estimate for the absolute pseudoranges at the previous epoch available, another piece of information for estimating the pseudoranges at the current epoch can be constructed by appropriately combining the absolute information at the previous epoch with the relative information in the interval from the previous to the current epoch. This kind of absolute information at the previous epoch is indeed available, however in the form of the coordinates plus clock error rather than the pseudoranges directly. Now it is clear that what we will do is merely a conversion, namely from the coordinates/clock-error to the pseudoranges. This conversion is nothing but an inverse pseudorange positioning problem, in fact an even simpler one than the direct one. Note that this conversion is not merely converting coordinates and clock error to pseudoranges, the associating covariance and cross covariance conversion is also performed. The pseudorange can be calculated for each satellite still visible at the previous epoch from the
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4937
estimated parameter vector as follows: j 2 j 2 j 2 j xk−1 − xˆk−1 + yk−1 pˆ k−1 = − yˆk−1 + zk−1 − zˆk−1 + tˆk−1
(26)
with the conversion of the corresponding covariance and cross covariance as follows: Q pˆ pˆ,k−1 = pˆcov pˆ k−1 = J T Qβˆ β,k−1 J ˆ
(27)
⎡
Hk−1 = cov pˆ k−1 , ϕk−1 = cov⎣J T βˆk−1 ,
⎤ j ⎦ τ j εk−1 = JT
j
j uk−1 τ Tj
(28)
j
where the column of the Jacobian J is aj , namely the coefficient vector corresponding to the jth satellite; and τ j’ denote an n × 1 vector with the j’th element being one and other elements being all zero, here n is the number of satellites still visible at the current epoch j and j’ denotes the index of pˆ k−1 in the vector pˆ k−1 , note that j’ may not necessarily be equal to j. If the carrier phase measurement from a satellite is involved only in the current epoch j but not in the previous epoch, the corresponding uk−1 in Eq. (28) is assigned zero. So the propagated pseudoranges from the previous epoch through time differenced carrier phases are as follows: p¯ k = pˆ k−1 + δϕk
(29)
T Q p¯ p¯,k = Q pˆ pˆ,k−1 + Qεε,k + Qεε,k−1 − Hk−1 − Hk−1
(30)
Finally we have the smoothed pseudoranges by simply performing the following leastsquares estimation:
−1 −1 Q pˆ pˆ,k = Q−1 (31) p¯ p¯,k + Qee,k
−1 ¯ p pˆ = Q pˆ pˆ,k Q−1 + Q ρ p¯ p¯,k k ee,k k
(32)
Hk = cov pˆ k , ϕk = cov Q pˆ pˆ,k Q−1 ϕ , ϕ = Q pˆ pˆ,k Q−1 k k p¯ p¯,k p¯ p¯,k Qε ε ,k
(33)
This completes the presentation of the transition algorithm from the PD to the global RD CSCF. 3.2. Global range-domain filter The global RD-CSCF at consecutive satellite-deficiency epochs is presented in the following. As stated before, all pseudoranges are processed as a whole to take care of the correlations, and hence the output is a smoothed pseudorange vector as at the PD-RD transition epoch stated above. At any one of consecutive deficient epochs, the input is exactly the output at the previous epoch, in the form as in Eqs. (31)–(33), namely the smoothed pseudorange vector and the corresponding covariance and cross covariance. First of all, it may be necessary to reconstruct the pseudorange vector and hence to reformulate the corresponding covariance and cross covariance matrix/vectors accordingly, due to potential visible satellite change. If a satellite becomes no longer visible, we simply delete the corresponding element
4938
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
from the vector as in Eq. (32) and the corresponding row and column in the matrix (31). Of course the column vector as in Eq. (33) corresponding to this satellite would be of no use and hence is deleted. If a new satellite becomes visible, however without ensuring the satellite sufficiency, we simply add the code pseudorange corresponding to this satellite to the final smoothed pseudorange vector as in Eq. (32), and assign the corresponding diagonal element as the variance of the differential code pseudorange. Due to the independence of this code pseudorange with other smoothed pseudoranges, all corresponding covariances are zero. This treatment of new visible satellite is the same as in the conventional RD-CSCF. Note that after incorporating the new visible satellite, the total number of the visible satellites is still less than four, otherwise, the algorithm for transition from satellite deficiency (global RD) to sufficiency (PD) will be triggered which is detailed in the next paragraph. For brevity, the variables after the reconstruction/reformulation are still denoted in the same form as the variables in Eqs. (31)–(33). The propagated/”predicted” pseudoranges and the corresponding covariance matrix are simply in the same form as in Eqs. (29) and (30). Then the final smoothed pseudoranges are obtained by fusing the propagated pseudoranges and the code pseudoranges in the same form as in Eqs. (31)–(33). However note that the dimensions and the corresponding satellites in the above expressions may be different from the previous epoch. 3.3. Switching from global range-domain to position-domain filters The transition algorithm from global RD at deficiency epoch to PD CSCF at a restoredsufficiency epoch is detailed. It is the coordinates/clock-error that is to be estimated at this transition epoch. The number of code pseudoranges is larger than or equal to four. However, the dimension of smoothed pseudoranges at the previous epoch is less than four, making any feasible coordinates/clock-error estimates unavailable. This fact shows that the conventional PD-CSCF algorithm presented in the previous section can no longer be applicable. The PDCSCF in the phase-connected-code form can still work [11]. However, for easy presentation and better linking with the output of the global RD algorithm at the previous deficiency epoch, the method in [11] is slightly modified. First, the propagated pseudoranges are obtained as in the same form as in Eqs. (29) and (30). Second, the smoothed pseudoranges is obtained by fusing the propagated pseudoranges and the corresponding code pseudoranges as in Eqs. (31) and (32). Note that the cross covariance matrix as in Eq. (33) is no longer needed. Third, stack the above smoothed pseudoranges with the additional code pseudoranges from other visible satellites as the observation vector, denoted as zk . The error covariance matrix, denoted as Qzz , k , for this observation vector is a block diagonal one, with the block corresponding to the smoothed pseudoranges being a full matrix calculated as in Eq. (31) and the other block corresponding to the additional code pseudorange/pseudoranges from other visible satellite/satellites being a scalar/diagonal matrix. Fourth, a least-squares pseudorange positioning is conducted in the similar for as in Eqs. (14) and (15) as follows: −1 T −1 βˆk = ATk Q−1 Ak Qzz,k zk = Sk zk (34) zz,k Ak T −1 −1 Qβˆ β,k ˆ = Ak Qzz,k Ak
(35)
2 ukj = cov βˆk , ϕkj = cov sk, j ϕkj , ϕkj = sk, j σϕ,k, j
(36)
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4939
where sk , j denotes the column vector of Sk corresponding to the jth satellite. It is not difficult to prove that the method presented here is equivalent to the carrier-connected-code method proposed in [11], though in different forms. Note that the cross covariance is only performed for the jth satellite whose propagated pseudorange, rather than merely the code pseudorange, is available; those for the satellites corresponding to which only the code pseudoranges are available are zero. The number of the time differenced carrier phases available at the current epoch may be larger than that of the smoothed pseudoranges at the previous epoch. However, only those whose corresponding smoothed pseudoranges at the previous epoch are available can be used here. Other time differenced carrier phases have to be abandoned at the current epoch. This is a defect of the proposed method owing to information loss; maybe future studies can be conducted to fix this problem. Before ending the presentation of the proposed switching PD and RD CSCF method, a flow chart is provided in Fig. 1 to show the detailed implementation algorithm of this method. Note that the case considered in Fig. 1 is more general than that dealt with in the above method presentation in that the agorithm may initializae at deficiency epoch. 4. Experiments Three algorithms are checked in the experiments, namely the conventional PD-CSCF, the conventional RD-CSCF and the proposed switching PD and RD CSCF, denoted SPD. A real BDS data set is analyzed. The static data are collected by ZHD receivers and V8 antenna in Huainan, Anhui Province. The baseline length is nearly 15 km, the sampling interval time is 1 s, and the epoch number is 1400. A cut-off elevation of 10° is set in all satellite data processing. All satellites of BDS, including GEO, IGSO, and MEO satellites, are used. The RTK results using all available measurements are employed as ground truth; and the errors of the east/north/up coordinates estimated in the three algorithms are calculated as the difference between them and the ground truth. The empirical variances for code pseudoranges and carrier phases from the zenith direction are set as 0.52 m2 and 0.0052 m2 , respectively. Without loss of generality, the variances along slant directions are obtained by inflating those along the zenith direction by the following factor: ⎧ π ⎪ for Ae ≥ ⎨1 , 6 ω= (37) π ⎪ ⎩2 sin (Ae ), for Ae < 6 where Ae denotes the elevation angle. The satellite deficiencies at some epochs are artificially introduced. According to the different kinds of deficiency, the following 8 different cases are checked. Case Case Case Case Case Case Case Case
1: 2: 3: 4: 5: 6: 7: 8:
3 2 1 1 3 2 1 1
visible visible visible visible visible visible visible visible
satellites in interval T1 ; satellites in interval T1 ; satellite in interval T1 ; satellite in intervals T1 and T2 ; satellites in interval T1 and 2 in T2 ; satellites in interval T1 and 1 in T2 ; satellite in interval T1 and 3 in T2 ; satellite in interval T1 , 2 in T2 and 3 in T3 .
4940
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
Visible satellites become insufficient
Visible satellites become sufficient
Fig. 1. Block diagram (Top) and flow chart (Bottom) of the switching PD and RD CSCF algorithm. PD, RD and CSCF denote position-domain, range-domain and carrier-smoothing-code filtering, respectively.
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4941
Fig. 2. Position error in East (Top), North (Middle) and Up (Bottom) directions for Case 1 (Left) and partial enlargement for satellite deficiencies epochs (Right), respectively.
All the above deficiency intervals are of 100-epoch length. The positioning errors are depicted in Figs. 2-4 for Case 1–3, respectively; these for case 8 are depicted in Figs. 5-7, for the three directions, respectively. The results are similar for other cases. Note that in these figures, the intervals, without any positioning solutions, are also depicted in the figures, intentionally to clearly show the complete time line. The detailed positioning errors of the PD and SPD algorithm at the 100 epochs immediately after a deficiency interval are displayed separately in the bottom part of all the 4 figures. The statistics of the errors for different cases is shown in Table 1. From all these figures and the table, the following can be observed: First of all, the superiority of the proposed SPD algorithm to conventional PD or RD algorithms is clearly demonstrated. To be more specific, both PD and SPD produce much smoother estimation than RD, this is because global optimality can be achieved by the former two and the latter, with a single-channel basis, can only be locally optimal. Although with similar steady state performance after a long period of continuous satellite sufficiency, SPD clearly outperforms PD with much improved transient accuracies at several epochs immediately after an interval of satellite deficiency. This is exactly what is expected by conducting this study. Here we stress again the mechanism of the superiority of SPD to PD: at epochs after an interval of satellite deficiency, historical information before and in this interval, completely abandoned in conventional PD, is effectively preserved and exploited in the proposed SPD. Another important observation can be made: estimation errors at the epoch immediately after a satellite-deficiency interval are much larger than that at the epoch immediately before this same interval. This is also reasonable as explained in the following: in any satellite deficiency interval, what we are doing with SPD is nothing but propagating the information in a
4942
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
Fig. 3. Position error in East (Top), North (Middle) and Up (Bottom) directions for Case 2 (Left) and partial enlargement for satellite deficiencies epochs (Right), respectively.
Fig. 4. Position error in East (Top), North (Middle) and Up (Bottom) directions for Case 3 (Left) and partial enlargement for satellite deficiencies epochs (Right), respectively.
NO.
1 2 3 4 5 6 7 8
E
N
U
T OTAL
Raw
RD
PD
SPD
Raw
RD
PD
SPD
Raw
RD
PD
SPD
Raw
RD
PD
SPD
0.530 0.530 0.530 0.530 0.534 0.534 0.534 0.536
0.190 0.190 0.190 0.190 0.199 0.199 0.199 0.196
0.133 0.133 0.133 0.170 0.197 0.197 0.197 0.186
0.118 0.121 0.123 0.133 0.168 0.171 0.170 0.154
0.855 0.855 0.855 0.855 0.890 0.890 0.890 0.893
0.315 0.315 0.315 0.315 0.324 0.324 0.324 0.321
0.275 0.275 0.275 0.298 0.318 0.318 0.318 0.305
0.234 0.239 0.242 0.275 0.276 0.282 0.282 0.265
0.452 0.451 0.451 0.499 0.499 0.499 0.499 0.499
0.195 0.195 0.195 0.220 0.220 0.220 0.220 0.199
0.190 0.190 0.190 0.216 0.216 0.216 0.216 0.219
0.178 0.184 0.185 0.209 0.201 0.207 0.206 0.206
1.126 1.102 1.102 1.123 1.131 1.152 1.152 1.155
0.416 0.416 0.416 0.429 0.439 0.439 0.439 0.426
0.360 0.360 0.360 0.402 0.432 0.432 0.432 0.419
0.317 0.325 0.329 0.374 0.381 0.389 0.388 0.369
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
Table 1 Mean squared errors for different cases (unit: meter).
4943
4944
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
Fig. 5. Position error in East directions for Case 8 (Top) and partial enlargement for satellite deficiencies epochs (Bottom), respectively.
Fig. 6. Position error in North directions for Case 8 (Top) and partial enlargement for satellite deficiencies epochs (Bottom), respectively.
four dimensional parameter vector estimate, namely the coordinates and clock error, however using a three, two, or one dimensional pseudorange estimate, namely the smoothed pseudorange(s). The dimension deficiency will inevitably cause information loss in the propagation. It can be expected that the longer the interval is, the more information there is, and hence the larger the errors after the interval than that before the interval. So with longer and longer satellite interval, the superiority of SPD to PD would be more and more marginal. For very long intervals, SPD would degenerate to PD. Seemingly, for this kind of case we prefer PD to
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4945
Fig. 7. Position error in Up directions for Case 8 (Top) and partial enlargement for satellite deficiencies epochs (Bottom), respectively.
SPD due to the simplicity of the former. However, we still recommend SPD for all cases with intermittent satellite deficiencies for the following two reasons: 1) it highly depends whether or not the interval is indeed long enough. It is apparently a difficult task to make such a decision, especially in real time and adaptively by taking all relevant factors into account. However, the prime aim of conducting this study is to deal with the case with short but intermittently occurred satellite deficiency intervals. Only for this kind of case, the performance improvement by employing the proposed SPD would be significant. 2) The complexity of the SPD at satellite deficiency epochs, namely the global RD, is low. In fact it is only a low-dimension recursive least squares problem with 3/2/1 parameters and 6/4/2 measurements. The only difference from an ordinary recursive least-squares estimation is the consideration of the special kind of coloredness of the measurement noises (the time differenced carrier phase measurements only); or to be more specific, only a cross covariance is additionally calculated, stored and update. A final note is that, for an interval without any visible satellites, none of the CSCF techniques can perform any information propagation and historical information has to be abandoned. This is reasonable, because for GNSS-only applications without introducing external aid, on media is available to carry the historical information. With this kind of interval, both conventional RD and the proposed SPD can re-initialize at an epoch when even one satellite becomes visible; however conventional PD has to re-initialize only at an epoch with restored satellite sufficiency. 5. Conclusion CSCF is an effective method to process the two basic types of GNSS signals: code pseudorange and carrier phase. This study aims to improve the positioning performance of GNSS-only CSCF technique under harsh environments with intermittent satellite deficiencies. By GNSSonly we mean that no aid is introduced, either from external sensors or from internal models. By satellite deficiency, we mean that there are visible satellites at some epochs; however the
4946
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
number is less than four. The PD filtering is superior to the RD one in terms of accuracy and insensitivity of satellite changes; however, the former is not applicable at satellite deficiency epochs. So the positioning performance improvement in this study is achieved by employing the generally superior PD filtering, however with special modifications tailored to consider the intermittently occurred satellite deficiencies. This modification is a seemingly simple switching, namely from the PD filtering for sufficiency epoch to the RD filtering for deficiency epoch and conversely. However, in the proposed switching method, the RD filtering for the deficiency epochs is a globally optimal one with all pseudoranges treated together as a vector in contrary to the conventional single-channel RD filtering which is only locally optimal. This global RD filter is necessary to preserve the correlations among different channels. This correlation indeed exist as at the first epoch of a deficiency interval, the smoothed pseudoranges is calculated using the same parameter vector estimate at the previous epoch. Real BDS data is processed and the superiority of the proposed method is validated, namely with higher steady state and transient accuracies than conventional RD and PD filtering, respectively. Acknowledgments This research is supported by the Fundamental Research Funds for the Central Universities of China under grant 2018QNA22. References [1] Z. Zhou, B. Li, Optimal Doppler-aided smoothing strategy for GNSS navigation, GPS Solut. 21 (2017) 197–210. [2] R.R. Hatch, The synergism of GPS code and carrier measurements, in: Third International Geodetic Symposium on Satellite Doppler Positioning, New Mexico, 1982, pp. 1213–1232. [3] P.J.G. Teunissen, The GPS phase-adjusted pseudorange, 2nd international workshop on high precision navigation, Stuttgart, 1991, pp. 115–125. [4] C. Gunther, P. Henkel, Reduced-noise ionosphere-free carrier smoothed code, IEEE Trans. Aerosp. Electron. Syst. 46 (2010) 323–334. [5] G. Chang, T. Xu, Y. Yao, Q. Wang, Adaptive Kalman filter based on variance component estimation for the prediction of ionospheric delay in aiding the cycle slip repair of GNSS triple-frequency signals, J. Geod. 92 (2018) 1241–1253. [6] P.Y.C. Hwang, G.A. McGraw, J.R. Bader, Enhanced differential GPS carrier-smoothed code processing using dual-frequency measurements, Navigation 46 (1999) 127–137. [7] G.A. McGraw, Generalized divergence-free carrier smoothing with applications to dual frequency differential GPS, Navigation 56 (2009) 115–122. [8] B. Park, C. Lim, Y. Yun, E. Kim, C. Kee, Optimal Divergence-free hatch filter for GNSS single-frequency measurement, Sensors 17 (2017) 448. [9] J.Y. Lee, H.S. Kim, K.H. Choi, J. Cho, H.K. Lee, Adaptive position-domain carrier-smoothed code filter based on innovation sequences, IET Radar, Sonar Navig. 8 (2013) 336–343. [10] H.K. Lee, C. Rizos, Position-domain hatch filter for kinematic differential GPS/GNSS, IEEE Trans. Aerosp. Electron. Syst. 44 (2008) 30–40. [11] S.B. Bisnath, R.B. Langley, Precise orbit determination of low earth orbiters with GPS point positioning, in: 2001 National Technical Meeting of the Institute of Navigation, Long Beach, CA, 2001, pp. 725–733. [12] H.K. Lee, C. Rizos, G.-I. Jee, Design of kinematic DGNSS filters with consistent error covariance information, IEE Proc.-Radar, Sonar Navig. 151 (2004) 382–388. [13] H.K. Lee, C. Rizos, G.-I. Jee, Position domain filtering and range domain filtering for carrier-smoothed-code DGNSS: an analytical comparison, IEE Proc.-Radar, Sonar Navig. 152 (2005) 271–276. [14] W. Liu, X. Wang, Z. Deng, Robust centralized and weighted measurement fusion Kalman estimators for multisensor systems with multiplicative and uncertain-covariance linearly correlated white noises, J. Frankl. Inst. 354 (2017) 1992–2031.
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4947
[15] L. Jiang, L. Yan, Y. Xia, Q. Guo, M. Fu, L. Li, Distributed fusion in wireless sensor networks based on a novel event-triggered strategy, J. Frankl. Inst. (2018) In Press. [16] G. Wang, N. Li, Y. Zhang, Maximum correntropy unscented Kalman and information filters for non-Gaussian measurement noise, J. Frankl. Inst. 354 (2017) 8659–8677. [17] W. Li, Y. Jia, J. Du, Distributed extended Kalman filter with nonlinear consensus estimate, J. Frankl. Inst. 354 (2017) 7983–7995.
Journal of the Franklin Institute 356 (2019) 4928–4947 www.elsevier.com/locate/jfranklin
Switching position and range-domain carrier-smoothing-code filtering for GNSS positioning in harsh environments with intermittent satellite deficiencies Guobin Chang a,c,∗, Tianhe Xu b,c, Chao Chen a, Bing Ji d, Shengquan Li e a School
of Environmental Science and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China b Institute of Space Science, Shandong University, Weihai 264209, China c State Key Laboratory of Geo-information Engineering, Xi’an Research Institute of Surveying and Mapping, Xi’an 710054, China d Department of Navigation Engineering, Naval University of Engineering, Wuhan 430033, China e Acoustic Science and Technology Laboratory, College of Underwater Acoustic Engineering, Harbin Engineering University, Harbin 150001, China Received 10 May 2018; received in revised form 19 December 2018; accepted 2 April 2019 Available online 11 April 2019
Abstract Carrier-smoothing-code filtering (CSCF) is widely used in GNSS signal processing to combine code pseudoranges and carrier phases. Position-domain (PD) CSCF is generally more accurate and less sensitive to visible satellite changes than range-domain (RD) CSCF. However, PD-CSCF necessitates at least four visible satellites. Intermittent satellite deficiency with less than four visible satellites is not uncommon in harsh environments like urban canyons. At such deficiency epochs, the PD-CSCF convergence has to break off. This study aims to bridge intermittent deficiency epochs in PD-CSCF without introducing any external information. The proposed solution is called switching RD and PD CSCF in which PD filter is replaced by RD filter at deficiency epochs. Besides detailing the seamless switching algorithms from PD/RD to RD/PD filters, a global RD filter, different from the conventional one, is developed to preserve the correlations among smoothed pseudoranges corresponding to different satellites. Compared to the conventional PD and RD CSCF algorithms, smoother results can be expected from the proposed switching filter, especially after the intermittent deficiency epochs. Experiments are ∗ Corresponding author at: School of Environmental Science and Spatial Informatics, China University of Mining and Technology. E-mail address: [email protected] (G. Chang).
https://doi.org/10.1016/j.jfranklin.2019.04.005 0016-0032/© 2019 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4929
conducted using real BDS signals. Cases with different kinds of deficiency are considered. Superiority of the proposed method is clearly observed from the results. © 2019 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.
1. Introduction Global Navigation Satellite Systems (GNSS) technology is widely used in broad engineering applications relating navigation, guidance and control. Basically, two types of signals are normally used in GNSS, namely code pseudoranges and carrier phases, though Doppler signals are also incorporated in some applications [1]. Compared with code pseudorange, carrier phase is with much higher accuracy but ambiguous with unknown integer cycles. In order to make full use of the high-accuracy carrier phase signals, two approaches can be used to address the ambiguity problem, namely performing integer ambiguity resolution (AR) including both estimation and validation and smoothing the code pseudoranges with carrier phases [2]. The latter, called the carrier-smoothing-code filtering (CSCF) technique, is our focus in this work. If we include integer cycles as parameters to be estimated as float values without fixing them as integers, the corresponding approach can roughly be regarded as a CSCF technique, though it is called carrier adjusted code (CAC) in [3]. In the CAC approach when applied in dynamical positioning, the integer cycles function as kind of media to carry historical information forward to current and future epochs. Of course, AR and CSCF approaches can be merged, e.g., before integer cycles can be fixed with sufficiently high success rate, CSCF can be used in the transient stage, or in cases of multiple-frequency signals, one can perform partial resolution of the integer ambiguities, namely linear combinations of the original integer ambiguities, with the remaining part dealt with using CSCF technique [4]. In CSCF, what one does is nothing but a combination of the code pseudoranges and the time differenced or epoch differenced carrier phases. The aim of the time difference operation is to cancel integer ambiguities as long as cycle slips are absent or repaired. Cycle slips detection and correction is a prerequisite for almost all applications involving carrier phases, not only for the CSCF studied in this work, but also for AR. Cycle slip detection/correction can and also have to be performed at different stages of the whole process of data processing, it is far from a trivial task, especially when performed in real time [5]. However, to be more focused, the cycle slip repairing issue is not to be explored in this study and we simple assume the absence of any cycle slips in the used carrier phase data. Time differenced carrier phases, carrying only relative or incremental information of the position, can be solely relied on in determining displacements, velocities or accelerations. However, they have to be combined with code pseudoranges to obtain non-divergent absolute positioning solutions, similarly as dead-reckon positioning techniques like inertial navigation system. The combination of code pseudoranges and time differenced carrier phases in the CSCF can be carried out generally in two different ways, namely in range domain (RD) and in position domain (PD). In RD-CSCF, signals are processed separately for different satellites [2,6–8]. However, in PD-CSCF, all available measurements, namely the code pseudoranges and the time differenced carrier phases from all visible satellites (considering the minimum/mask elevation angle), are processed altogether to estimate the position coordinates and of course the receiver clock error. This is nothing but exactly a process of information fusion. Although other fusion strategies or some extra tunings
4930
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
may be employed [9,10], the fusion is often optimally performed in the least-squares sense [11–13]. When used in static positioning, e.g., at the reference stations of ground or satellite based augmentation systems, this filtering degenerates to a recursive least-squares estimation. It has been shown, both theoretically and experimentally, PD-CSCF can be superior to RDCSCF in general in terms of accuracy and insensitivity to visible satellite changes [10,13]. This is not any surprise when considered in the framework of information fusion or statistical signal processing. In PD-CSCF, besides the information propagation along time line, further information sharing across different satellites/channels is also permitted. Satellite redundancy, not merely sufficiency, namely with more than 4 visible satellites, can be exploited fully by PD-CSCF which RD-CSCF fails to do. The more visible satellites there are, the more obvious of the advantage of the PD-CSCF will be. Note that the following needs to be emphasized to avoid misunderstanding: with exactly 4 satellites visible, namely without redundancy, PDCSCF and RD-CSCF would not be equivalent to each other either, in general, or to be more specific, the former would still outperform the latter. The more accurate PD-CSCF is also more restrictive. At least four satellites have to be visible in order to perform the PD-CSCF. As the single-channel RD-CSCF is carried out separately for each visible satellite, there is no such restriction on the number of visible satellites, though positioning using the output of the RD-CSCF, namely the smoothed pseudoranges, can be possible only with at least four visible satellites. For scenarios with intermittent satellite deficiencies, RD-CSCF would outperform PD-CSCF in a sense as explained in the following. Although neither can provide positioning service at the deficiency epochs, RD-CSCF can still continuously pass the historical information buried in the still visible satellites on to the smoothed pseudoranges at these deficiency epochs and also on to those at future sufficiency epochs; however, PD-CSCF can NOT ensure such continuous information passing. Satellite deficiency at any epoch would break off the converging process of the PD-CSCF. The historical information before the deficiency epoch has to be entirely abandoned, and the convergence has to start over again when satellite sufficiency restores. Frequent deficiencies would break the filtering process into pieces in which smoothing effects could not be obvious owing to the short converging periods. This of course would result in poor positioning performances. Scenarios of satellite deficiency are not uncommon in some harsh environments with blocking/jamming or multipath, such as urban canyons. In the past, this deficiency is mitigated mainly by fusing additional information [14,15]. The additional information can be those provided by external sensors and/or infrastructures. Exploiting the model of the environment, the vehicle or the equipment/sensor, such as the road maps, the vehicle kinematics models, and the receiver clock models can also provide information to permit a solution during satellite deficiencies. Of course these efforts are meaningful, however at a price, e.g., cost and/or complexity would be increased. There are some situations in which short period of positioning unavailability can be tolerated, e.g., when a car drives under a viaduct or passes by a tall building. In these situations, cumbersome job of adding external sensors or exploiting models may not be cost effective. In this work it is tried to generalize PD-CSCF to the case with intermittent satellite deficiencies and without any external information. The solution, seemingly straightforward, can be using a switching PD and RD CSCF approach, namely using PD-CSCF for satellitesufficiency epochs and RD-CSCF for satellite-deficiency epochs. Then it is worthwhile to conduct a seamless transition from PD-CSCF to RD-CSCF and of course also a seamless opposite transition. By saying seamless, we mean that with the transition from PD to RD, all relevant statistical information in the coordinates/clock-error estimate in the PD-CSCF would
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4931
be propagated into the smoothed pseudoranges in the RD-CSCF and vice versa. To our best knowledge, this switching strategy has not been studied in the literature. Our aim to conduct this study is to present the detailed model and algorithm for the switching PD and RD CSCF for GNSS-only applications in cases with intermittent satellite deficiencies. While the PD algorithm for epochs of satellite sufficiency in the proposed switching PD and RD CSCF method is exactly the same to the conventional PD-CSCF, the RD algorithm for epochs of satellite deficiency deviates from the traditional RD-CSCF algorithm. At the transition epoch from satellite sufficiency to deficiency, the same set of statistical information at previous epoch flows into any one of the smoothed pseudoranges at the current epoch. Apparently, the thus obtained smoothed pseudoranges are no longer independent in contrary to the conventional channel-wise RD-CSCF. At the subsequent epochs with satellite deficiency, if we conduct the conventional single-channel RD filtering for each visible satellite separately as in the conventional RD-CSCF, the solution would not be statistically optimal owing to ignoring the correlations among different channels. So instead, in the proposed switching algorithm, the RD processing at deficiency epochs for all visible satellites are performed as a whole rather than separately. This is statistically correct and brings in global optimality. At deficiency epochs, both the switching method and the conventional RD-only method are conducting RD filtering, however, each one of the smoothed pseudoranges provided by the latter only contain the historical information of the corresponding satellite while the former also contain the historical information of other, previously visible satellites. This would produce smoother results enabling higher accuracy and faster convergence. This is exactly what we are pursuing by conducting this study. The single-frequency point positioning is employed as an example in both the method presentation and the experiment study. Assume all error terms except clock error and noises have been corrected using models. This simplifies the model with most error terms absent and makes us more focused on the key element, namely the switching PD and RD CSCF algorithm development and validation. This is of course without loss of generality; the proposed method can be readily extended to other cases, e.g., relative (differential) positioning and/or multiple frequencies. For relative positioning with short to mediate baselines, the problem becomes even easier with most common error terms being cancelled or reduced largely [8]; and for multiple frequencies, linear combinations of raw measurements may be used as the inputs of the CSCF techniques including the one proposed in this work [4,6,7]. Of course, the combination and the CSCF can be performed in one step where global optimality can be achieved by selecting the combining and filtering coefficients simultaneously. The rest of the paper is organized as follows. In Section 2, the conventional PD and RD CSCF algorithms are presented in order to introduce terminology and to make the proposed method in the context. The main algorithm is presented in Section 3. Experiments using real BDS data are presented in Section 4 and the results are analyzed to check the performance of the proposed method. Concluding remarks are given in section V also with some potential topics of future study suggested.
2. Range-domain and position domain carrier-smoothing-code filters Kinematic point positioning is employed here. Error terms are assumed to be absent or removed/corrected using empirical models, such as satellite orbit and clock errors, tropospheric and ionospheric delay errors.
4932
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
2.1. Measurement model The code pseudorange and carrier phase (in meters) measurement equation would then be in the following simple forms [3]: ρkj = rkj + tk + ekj
(1)
ϕkj = rkj + tk + λnk + εkj
(2)
In the above, subscript k and superscript j denote the indices of epoch and the satellite, respectively; ρ and ϕ denote the code pseudoranges and carrier phases, respectively and e and ε denote the corresponding noises including all un-modeled errors, such as hardware noises and multipath errors, satellite orbit and clock errors and residual tropospheric/ionospheric delay errors; λ and n denote the carrier frequency and the unknown integer cycles, respectively; t denotes the receiver clock error multiplied by the light speed; r denotes the geometric range whose expression is as follows: j 2 2 2 xk − xk + ykj − yk + zkj − zk rkj = (3) where x/y/z denotes the Cartesian coordinate triad; the x/y/z with superscript j denote the known coordinates of the satellite and those without any superscript denote the coordinates of the rover receiver which, together with tk , are the parameters to be estimated. In most CSCF approaches, time differencing carrier phases is firstly conducted to remove the unknown but constant integer cycles; so the following measurement equation is used instead of Eq. (2): j δϕkj = δrkj + δtk + εkj − εk−1
(4)
Of course Eq. (4) relies on the absence of cycle slips. Preprocessing of detecting and repairing the potential cycle slips is needed, which is out of the range of this study. 2.2. Range-domain carrier-smoothing-code filters In RD-CSCF, the range plus receiver clock error, namely, pkj = rkj + tk
(5)
called (true) pseudorange, is the parameter firstly to be estimated, though the coordinates of the rover station are of final interest. The first and most widely used version of this kind of approach is the Hatch filter, which is a first-order infinite impulse response filter with fixed or varying but predetermined filtering coefficients [2]. However, in this work, a stochastic filtering methodology based on the least-squares criterion, rather than the digital filtering in the original Hatch filter, is followed. There are no essential differences between the stochastic and the Hatch filters; we are only replacing the predetermined filtering coefficient in the latter with a variant one adaptively tuned in the former. This variant coefficient is dynamically calculated according to the variances of the involved variables. The filtering algorithm is presented for the initialization and the recursion stages separately for convenience. For the initialization stage, the output, namely the smoothed pseudorange, as kind of estimate of the (true) pseudorange, denoted as pˆ kj , is assigned directly as the code pseudorange. Of course the variances of the smoothed pseudoranges are assigned as those of the code pseudoranges. Then for the case at the initialization epoch, the covariance defined in the following is zero: ckj = cov pˆ kj , ϕkj = 0
(6)
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4933
This is clearly reasonable, because there is not any carrier phase signals involved in calculating pˆ kj at the initialization epoch. Note that the epoch k does not necessarily equate to zero or one. Any epoch at which a satellite starts to be visible is an initialization epoch for this satellite. Any epoch at which cycle slip is detected but failed to be corrected (considering the success rate) is also viewed as an initialization epoch. The algorithm for the recursion is given in the following. For each epoch, there are two kinds of measurement for estimating the pseudorange. The first is that propagated from the previous epoch using the time differenced carrier phases, namely: j j j p¯ kj = pˆ k−1 +δϕkj = pkj + ξk−1 + εkj − εk−1 = pkj + ςkj j ξk−1
(7)
j pˆ k−1
2 where denotes the error in with the corresponding variance denoted as σˆ k−1 , j . The variance of the error in the above measurement is as follows: j j j j + var εkj + var εk−1 − 2cov ξk−1 σ¯ k,2 j = var ςkj = var ξk−1 , εk−1 j 2 2 2 = σˆ k−1 , j + σϕ,k, j + σϕ,k−1, j − 2ck−1
(8)
j 2 2 2 where σˆ k−1 , j and ck−1 are calculated at the previous epoch and σϕ,k, j /σϕ,k−1, j are the variances of the carrier phase measurements which are assigned empirically and based on elevation angle. The second kind of measurement is the raw code pseudorange as shown in Eq. (1). Combine the measurements in Eqs. (1) and (7) optimally in the least-squares sense as follows:
j −2 pˆ kj = σˆ k,2 j σ¯ k,−2j p¯ kj +σρ,k, j ρk
(9)
with
−2 −1 σˆ k,2 j = var pˆ kj = σ¯ k,−2j +σρ,k, j
(10)
2 ck = cov pˆ kj , ϕkj = cov σˆ k,2 j σ¯ k,−2j ϕkj , ϕkj = σˆ k,2 j σ¯ k,−2j σϕ,k, j
(11)
2 where σρ,k, j is the variance of the code pseudorange measurements which is also assigned empirically and based on elevation angle. Note that besides the variances in Eq. (10), the cross covariance calculated in Eq. (11) is necessary for the next epoch as is apparent in Eq. (8). At any epoch, with the smoothed pseudoranges available from at least four satellites, the coordinates of the rover station can be determined together with the receiver clock error, in the same way as the standard code pseudorange positioning. Note that all error terms common to code pseudoranges and carrier phases, e.g., the satellite orbit and clock errors and the residual tropospheric delay errors, do not matter in the determination of the smoothed pseudoranges, because all these errors are absorbed in the pseudoranges. It is the coordinate determination in which these errors enter as factors; however, they are assumed small after the model correction in this work.
2.3. Position-domain carrier-smoothing-code filters Contrary to the two-step strategy in the RD-CSCF algorithm presented in the above, namely first smoothing pseudoranges and then determining coordinates/clock-error, the PD-CSCF employs a one-step strategy, namely, determining the coordinates/clock-error directly using all available code pseudoranges and time differenced carrier phases. Let βk = [xk yk zk tk ]T be the parameter vector. Similarly to the RD-CSCF, at the initialization epoch, the parameter vector
4934
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
is estimated solely using code pseudoranges, and the cross covariance matrix between the parameter estimates and the carrier phase measurement from the jth satellite at the same initialization epoch, denoted as ukj = cov[βˆk ,ϕkj ], is a 4 × 1 zero vector. Note that visible satellite set changes will not trigger re-initialization as long as there are still at least four satellites visible. This is exactly the characteristic that PD-CSCF is not sensitive to satellite changes as stated above. However, at any epoch with satellite sufficiency restoring from previous deficiency, the PD-CSCF will re-initialize and hence this epoch is also an initialization epoch for the conventional PD-CSCF. For detailed description of PD-CSCF see e.g., [10,13]. Substitute Eq. (3) into Eq. (1); perform linearization at some approximate values of the coordinate; and we have the following measurement equation: lkj = aTk, j βk + ekj
(12)
where l and a denote the observation and design coefficient, respectively, with the superscript and subscript defined as before. Note that in Eq. (12), we directly treat the parameter vector βk as the argument rather than the correction or delta vector, for convenient presentation of the method. Either the coordinates from the standard positioning using only code pseudoranges or from the solution at the previous epoch can be used as the approximate values for the linearization, namely calculating the a vectors in Eq. (12). Stack equations as in Eq. (12) for all visible satellites, we have the vector measurement equation: lk = Ak βk + ek
(13)
Then the least-squares estimation is as follows: −1 T −1 β˜k = ATk Q−1 Ak Qee,k lk ee,k Ak
(14)
T −1 −1 Qβ˜ β,k ˜ = Ak Qee,k Ak
(15)
This represents the first piece of information for the parameter. The second piece is derived as follows. Linearize Eq. (4) with respect to time differenced vector δβk as follows: j j T ykj = ak, j δβk + εk − εk−1
(16)
where y denote the observation; and a is the same as defined in Eq. (12). Stack all available equations as in Eq. (16), we have the vector measurement equation as follows: yk = Bk δβk + εk − εk−1
(17)
Note here though the coefficient a in Eq. (16) is the same as that in Eq. (12), the coefficient matrix B in Eq. (17) does not necessarily be the same as that in Eq. (13). The visible satellite set represented in Eq. (16) may not be the same as that in Eq. (12), because a satellite represented in the former has to be continuously tracked at both the previous and the current epochs, while that in the latter only needs to be tracked at the current epoch. The least-squares estimation for Eq. (17) is as follows: −1 δ β¯k = BTk Qεε,k + Qεε,k−1 Bk −1 BTk Qεε,k + Qεε,k−1 −1 yk = Fk yk (18) T −1 −1 Qδβδ Bk ¯ β,k ¯ = Bk Qε ε ,k + Qε ε ,k−1
(19)
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
Let fk,j denote the jth column of Fk . Then we have the following cross covariance: ⎡ ⎤ j j ⎦ Gk = cov βˆk−1 , δ β¯k = cov⎣βˆk−1 , − fk, j εk−1 =− uk−1 fk,T j j∈ 0
4935
(20)
j∈ 0
In the above 0 denotes the index set of all the visible satellites from which the time differenced carrier phase signals are both available at the (k − 1)th and the kth epochs. Or in other words, a satellite whose index is in 0 has to be continuously tracked at the (k − 2)th, the (k − 1)th and the kth epochs. Then we have the propagated estimates for the parameter at the kth epoch as follows: β¯k = δ β¯k + βˆk−1
(21)
Qβ¯ β,k + GkT + Gk ¯ = Qδ βδ ¯ β,k ¯ + Qβˆ β,k−1 ˆ
(22)
This represents the second piece of information for estimating the parameter at the current epoch. It is clear that these two pieces of information are independent. So by simply fusing them optimally in the least-squares sense, we have the following expression for the final parameter estimates at the current epoch:
−1 −1 −1 Q Qβˆ β,k = + Q (23) ˆ ¯ ˜ β¯ β,k β˜ β,k
−1 ¯ Qβ−1 βˆk = Qβˆ β,k β k ˆ ¯ β,k ¯ βk + Qβ˜ β,k ˜
(24)
j j −1 −1 2 ukj = cov βˆk , ϕkj = cov Qβˆ β,k Q f ϕ , ϕ k, j ˆ ˆ Qβ¯ β,k ¯ ¯ f k, j σϕ,k, j k k = Qβˆ β,k β¯ β,k
(25)
Again, similar to the RD-CSCF, besides the covariance matrix in (23), the cross covariance (vector) in (25) is also necessary which will be used in the estimation at the next epoch. This cross covariances exactly represents the specific consideration of the time correlations in the time difference carrier phase signals. The PD-CSCF algorithm presented here is exactly the same to that derived in [12,13], though with a slightly different presentation. The following two are noted concerning the PD-CSCF. First, the second piece of information, namely the estimate in Eq. (21) may be called a prediction in some references. However, this cannot be a true prediction as that in the conventional Kalman filter. As is clear from the above derivation, this estimate, namely the so called prediction, can only be available at the current epoch, the same to the estimate in Eq. (14), because the carrier phase signals at the current epoch are involved in constructing the time differenced carrier phase signals. This is contrary to the conventional Kalman filter in which the prediction or time update is carried out at the (k − 1)th epoch [16,17]. So in this sense, there is no prediction of any physical meaning at all in the PD-CSCF. Second, from the presentation given here, both the number of the code pseudoranges and that of the time differenced carrier phases have to be at least four. However, this is an excessive demand. In fact either of the two numbers being larger than or equal to four can be sufficient for the existence of a parameter estimate at the current epoch. Or even less stringently, the PD-CSCF algorithm will be applicable as long as measurements from at least four different satellites are available; no matter the measurements are the code pseudoranges or the time differenced carrier phases. This is clearer when the PD-CSCF algorithm is presented in the phase-connected-code form [11]. Finally after the presentation, we want to stress once
4936
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
more that the PD-CSCF, no matter in what form, necessitates that at least four satellites have to be visible above the predetermined minimum/mask elevation angle. The CSCF algorithm for cases with intermittent satellite deficiencies, which is the main contribution of this work, is developed in detail in the next section.
3. Switching range and position domain carrier-smoothing-code filter In this section, the switching PD and RD CSCF algorithm is presented. As stated before, the RD algorithm in the proposed switching method at deficiency epochs deviates from the conventional single-channel RD algorithm presented in the previous section; so the RD algorithm here, essentially a global one, is also detailed in this section. Without loss of generality, the switching algorithm is initialized for a satellite-sufficiency epoch as the conventional PDCSCF presented in the last part of the previous section. As a result, in the following, the switching algorithm is presented in three parts, namely for the transition from the PD to the global RD filters, for the global RD filtering at the intermittent satellite-deficiency epochs, and for the transition from the global RD to the PD filters, respectively.
3.1. Switching from position-domain to global range-domain filters Once the number of visible satellites becomes deficient, namely with sufficient visible satellites at the immediately previous epoch, the transition from the PD to the global RD CSCF is triggered. The input to this transition is exactly the same to normal PD algorithm, namely the parameter estimate, the associate covariance matrix and the cross covariance vectors between the parameter estimate and the involved carrier phases. At this transition epoch, the parameters to be estimated changed from the coordinates/clock-error to the pseudoranges. Stack all the pseudoranges corresponding to the still visible satellites as the pseudorange vector, denoted as pk . Note that for 2 or 3 visible satellites, this vector are indeed a vector, with 2 and 3 elements, respectively. However, for the case of only one visible satellite, this vector degenerates to a scalar. Similarly stack code pseudorange and carrier phase measurements for all satellite still visible at the current epoch, in the same order as in pk , and denote them as the code pseudorange vector ρk and the carrier phase vector ϕk , respectively. For estimating the pseudorange vector, the code pseudoranges provide a piece of information. The time differenced carrier phases provide a piece of information for the relative or incremental pseudoranges in the interval from the previous to the current epoch. With the estimate for the absolute pseudoranges at the previous epoch available, another piece of information for estimating the pseudoranges at the current epoch can be constructed by appropriately combining the absolute information at the previous epoch with the relative information in the interval from the previous to the current epoch. This kind of absolute information at the previous epoch is indeed available, however in the form of the coordinates plus clock error rather than the pseudoranges directly. Now it is clear that what we will do is merely a conversion, namely from the coordinates/clock-error to the pseudoranges. This conversion is nothing but an inverse pseudorange positioning problem, in fact an even simpler one than the direct one. Note that this conversion is not merely converting coordinates and clock error to pseudoranges, the associating covariance and cross covariance conversion is also performed. The pseudorange can be calculated for each satellite still visible at the previous epoch from the
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4937
estimated parameter vector as follows: j 2 j 2 j 2 j xk−1 − xˆk−1 + yk−1 pˆ k−1 = − yˆk−1 + zk−1 − zˆk−1 + tˆk−1
(26)
with the conversion of the corresponding covariance and cross covariance as follows: Q pˆ pˆ,k−1 = pˆcov pˆ k−1 = J T Qβˆ β,k−1 J ˆ
(27)
⎡
Hk−1 = cov pˆ k−1 , ϕk−1 = cov⎣J T βˆk−1 ,
⎤ j ⎦ τ j εk−1 = JT
j
j uk−1 τ Tj
(28)
j
where the column of the Jacobian J is aj , namely the coefficient vector corresponding to the jth satellite; and τ j’ denote an n × 1 vector with the j’th element being one and other elements being all zero, here n is the number of satellites still visible at the current epoch j and j’ denotes the index of pˆ k−1 in the vector pˆ k−1 , note that j’ may not necessarily be equal to j. If the carrier phase measurement from a satellite is involved only in the current epoch j but not in the previous epoch, the corresponding uk−1 in Eq. (28) is assigned zero. So the propagated pseudoranges from the previous epoch through time differenced carrier phases are as follows: p¯ k = pˆ k−1 + δϕk
(29)
T Q p¯ p¯,k = Q pˆ pˆ,k−1 + Qεε,k + Qεε,k−1 − Hk−1 − Hk−1
(30)
Finally we have the smoothed pseudoranges by simply performing the following leastsquares estimation:
−1 −1 Q pˆ pˆ,k = Q−1 (31) p¯ p¯,k + Qee,k
−1 ¯ p pˆ = Q pˆ pˆ,k Q−1 + Q ρ p¯ p¯,k k ee,k k
(32)
Hk = cov pˆ k , ϕk = cov Q pˆ pˆ,k Q−1 ϕ , ϕ = Q pˆ pˆ,k Q−1 k k p¯ p¯,k p¯ p¯,k Qε ε ,k
(33)
This completes the presentation of the transition algorithm from the PD to the global RD CSCF. 3.2. Global range-domain filter The global RD-CSCF at consecutive satellite-deficiency epochs is presented in the following. As stated before, all pseudoranges are processed as a whole to take care of the correlations, and hence the output is a smoothed pseudorange vector as at the PD-RD transition epoch stated above. At any one of consecutive deficient epochs, the input is exactly the output at the previous epoch, in the form as in Eqs. (31)–(33), namely the smoothed pseudorange vector and the corresponding covariance and cross covariance. First of all, it may be necessary to reconstruct the pseudorange vector and hence to reformulate the corresponding covariance and cross covariance matrix/vectors accordingly, due to potential visible satellite change. If a satellite becomes no longer visible, we simply delete the corresponding element
4938
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
from the vector as in Eq. (32) and the corresponding row and column in the matrix (31). Of course the column vector as in Eq. (33) corresponding to this satellite would be of no use and hence is deleted. If a new satellite becomes visible, however without ensuring the satellite sufficiency, we simply add the code pseudorange corresponding to this satellite to the final smoothed pseudorange vector as in Eq. (32), and assign the corresponding diagonal element as the variance of the differential code pseudorange. Due to the independence of this code pseudorange with other smoothed pseudoranges, all corresponding covariances are zero. This treatment of new visible satellite is the same as in the conventional RD-CSCF. Note that after incorporating the new visible satellite, the total number of the visible satellites is still less than four, otherwise, the algorithm for transition from satellite deficiency (global RD) to sufficiency (PD) will be triggered which is detailed in the next paragraph. For brevity, the variables after the reconstruction/reformulation are still denoted in the same form as the variables in Eqs. (31)–(33). The propagated/”predicted” pseudoranges and the corresponding covariance matrix are simply in the same form as in Eqs. (29) and (30). Then the final smoothed pseudoranges are obtained by fusing the propagated pseudoranges and the code pseudoranges in the same form as in Eqs. (31)–(33). However note that the dimensions and the corresponding satellites in the above expressions may be different from the previous epoch. 3.3. Switching from global range-domain to position-domain filters The transition algorithm from global RD at deficiency epoch to PD CSCF at a restoredsufficiency epoch is detailed. It is the coordinates/clock-error that is to be estimated at this transition epoch. The number of code pseudoranges is larger than or equal to four. However, the dimension of smoothed pseudoranges at the previous epoch is less than four, making any feasible coordinates/clock-error estimates unavailable. This fact shows that the conventional PD-CSCF algorithm presented in the previous section can no longer be applicable. The PDCSCF in the phase-connected-code form can still work [11]. However, for easy presentation and better linking with the output of the global RD algorithm at the previous deficiency epoch, the method in [11] is slightly modified. First, the propagated pseudoranges are obtained as in the same form as in Eqs. (29) and (30). Second, the smoothed pseudoranges is obtained by fusing the propagated pseudoranges and the corresponding code pseudoranges as in Eqs. (31) and (32). Note that the cross covariance matrix as in Eq. (33) is no longer needed. Third, stack the above smoothed pseudoranges with the additional code pseudoranges from other visible satellites as the observation vector, denoted as zk . The error covariance matrix, denoted as Qzz , k , for this observation vector is a block diagonal one, with the block corresponding to the smoothed pseudoranges being a full matrix calculated as in Eq. (31) and the other block corresponding to the additional code pseudorange/pseudoranges from other visible satellite/satellites being a scalar/diagonal matrix. Fourth, a least-squares pseudorange positioning is conducted in the similar for as in Eqs. (14) and (15) as follows: −1 T −1 βˆk = ATk Q−1 Ak Qzz,k zk = Sk zk (34) zz,k Ak T −1 −1 Qβˆ β,k ˆ = Ak Qzz,k Ak
(35)
2 ukj = cov βˆk , ϕkj = cov sk, j ϕkj , ϕkj = sk, j σϕ,k, j
(36)
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4939
where sk , j denotes the column vector of Sk corresponding to the jth satellite. It is not difficult to prove that the method presented here is equivalent to the carrier-connected-code method proposed in [11], though in different forms. Note that the cross covariance is only performed for the jth satellite whose propagated pseudorange, rather than merely the code pseudorange, is available; those for the satellites corresponding to which only the code pseudoranges are available are zero. The number of the time differenced carrier phases available at the current epoch may be larger than that of the smoothed pseudoranges at the previous epoch. However, only those whose corresponding smoothed pseudoranges at the previous epoch are available can be used here. Other time differenced carrier phases have to be abandoned at the current epoch. This is a defect of the proposed method owing to information loss; maybe future studies can be conducted to fix this problem. Before ending the presentation of the proposed switching PD and RD CSCF method, a flow chart is provided in Fig. 1 to show the detailed implementation algorithm of this method. Note that the case considered in Fig. 1 is more general than that dealt with in the above method presentation in that the agorithm may initializae at deficiency epoch. 4. Experiments Three algorithms are checked in the experiments, namely the conventional PD-CSCF, the conventional RD-CSCF and the proposed switching PD and RD CSCF, denoted SPD. A real BDS data set is analyzed. The static data are collected by ZHD receivers and V8 antenna in Huainan, Anhui Province. The baseline length is nearly 15 km, the sampling interval time is 1 s, and the epoch number is 1400. A cut-off elevation of 10° is set in all satellite data processing. All satellites of BDS, including GEO, IGSO, and MEO satellites, are used. The RTK results using all available measurements are employed as ground truth; and the errors of the east/north/up coordinates estimated in the three algorithms are calculated as the difference between them and the ground truth. The empirical variances for code pseudoranges and carrier phases from the zenith direction are set as 0.52 m2 and 0.0052 m2 , respectively. Without loss of generality, the variances along slant directions are obtained by inflating those along the zenith direction by the following factor: ⎧ π ⎪ for Ae ≥ ⎨1 , 6 ω= (37) π ⎪ ⎩2 sin (Ae ), for Ae < 6 where Ae denotes the elevation angle. The satellite deficiencies at some epochs are artificially introduced. According to the different kinds of deficiency, the following 8 different cases are checked. Case Case Case Case Case Case Case Case
1: 2: 3: 4: 5: 6: 7: 8:
3 2 1 1 3 2 1 1
visible visible visible visible visible visible visible visible
satellites in interval T1 ; satellites in interval T1 ; satellite in interval T1 ; satellite in intervals T1 and T2 ; satellites in interval T1 and 2 in T2 ; satellites in interval T1 and 1 in T2 ; satellite in interval T1 and 3 in T2 ; satellite in interval T1 , 2 in T2 and 3 in T3 .
4940
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
Visible satellites become insufficient
Visible satellites become sufficient
Fig. 1. Block diagram (Top) and flow chart (Bottom) of the switching PD and RD CSCF algorithm. PD, RD and CSCF denote position-domain, range-domain and carrier-smoothing-code filtering, respectively.
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4941
Fig. 2. Position error in East (Top), North (Middle) and Up (Bottom) directions for Case 1 (Left) and partial enlargement for satellite deficiencies epochs (Right), respectively.
All the above deficiency intervals are of 100-epoch length. The positioning errors are depicted in Figs. 2-4 for Case 1–3, respectively; these for case 8 are depicted in Figs. 5-7, for the three directions, respectively. The results are similar for other cases. Note that in these figures, the intervals, without any positioning solutions, are also depicted in the figures, intentionally to clearly show the complete time line. The detailed positioning errors of the PD and SPD algorithm at the 100 epochs immediately after a deficiency interval are displayed separately in the bottom part of all the 4 figures. The statistics of the errors for different cases is shown in Table 1. From all these figures and the table, the following can be observed: First of all, the superiority of the proposed SPD algorithm to conventional PD or RD algorithms is clearly demonstrated. To be more specific, both PD and SPD produce much smoother estimation than RD, this is because global optimality can be achieved by the former two and the latter, with a single-channel basis, can only be locally optimal. Although with similar steady state performance after a long period of continuous satellite sufficiency, SPD clearly outperforms PD with much improved transient accuracies at several epochs immediately after an interval of satellite deficiency. This is exactly what is expected by conducting this study. Here we stress again the mechanism of the superiority of SPD to PD: at epochs after an interval of satellite deficiency, historical information before and in this interval, completely abandoned in conventional PD, is effectively preserved and exploited in the proposed SPD. Another important observation can be made: estimation errors at the epoch immediately after a satellite-deficiency interval are much larger than that at the epoch immediately before this same interval. This is also reasonable as explained in the following: in any satellite deficiency interval, what we are doing with SPD is nothing but propagating the information in a
4942
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
Fig. 3. Position error in East (Top), North (Middle) and Up (Bottom) directions for Case 2 (Left) and partial enlargement for satellite deficiencies epochs (Right), respectively.
Fig. 4. Position error in East (Top), North (Middle) and Up (Bottom) directions for Case 3 (Left) and partial enlargement for satellite deficiencies epochs (Right), respectively.
NO.
1 2 3 4 5 6 7 8
E
N
U
T OTAL
Raw
RD
PD
SPD
Raw
RD
PD
SPD
Raw
RD
PD
SPD
Raw
RD
PD
SPD
0.530 0.530 0.530 0.530 0.534 0.534 0.534 0.536
0.190 0.190 0.190 0.190 0.199 0.199 0.199 0.196
0.133 0.133 0.133 0.170 0.197 0.197 0.197 0.186
0.118 0.121 0.123 0.133 0.168 0.171 0.170 0.154
0.855 0.855 0.855 0.855 0.890 0.890 0.890 0.893
0.315 0.315 0.315 0.315 0.324 0.324 0.324 0.321
0.275 0.275 0.275 0.298 0.318 0.318 0.318 0.305
0.234 0.239 0.242 0.275 0.276 0.282 0.282 0.265
0.452 0.451 0.451 0.499 0.499 0.499 0.499 0.499
0.195 0.195 0.195 0.220 0.220 0.220 0.220 0.199
0.190 0.190 0.190 0.216 0.216 0.216 0.216 0.219
0.178 0.184 0.185 0.209 0.201 0.207 0.206 0.206
1.126 1.102 1.102 1.123 1.131 1.152 1.152 1.155
0.416 0.416 0.416 0.429 0.439 0.439 0.439 0.426
0.360 0.360 0.360 0.402 0.432 0.432 0.432 0.419
0.317 0.325 0.329 0.374 0.381 0.389 0.388 0.369
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
Table 1 Mean squared errors for different cases (unit: meter).
4943
4944
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
Fig. 5. Position error in East directions for Case 8 (Top) and partial enlargement for satellite deficiencies epochs (Bottom), respectively.
Fig. 6. Position error in North directions for Case 8 (Top) and partial enlargement for satellite deficiencies epochs (Bottom), respectively.
four dimensional parameter vector estimate, namely the coordinates and clock error, however using a three, two, or one dimensional pseudorange estimate, namely the smoothed pseudorange(s). The dimension deficiency will inevitably cause information loss in the propagation. It can be expected that the longer the interval is, the more information there is, and hence the larger the errors after the interval than that before the interval. So with longer and longer satellite interval, the superiority of SPD to PD would be more and more marginal. For very long intervals, SPD would degenerate to PD. Seemingly, for this kind of case we prefer PD to
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4945
Fig. 7. Position error in Up directions for Case 8 (Top) and partial enlargement for satellite deficiencies epochs (Bottom), respectively.
SPD due to the simplicity of the former. However, we still recommend SPD for all cases with intermittent satellite deficiencies for the following two reasons: 1) it highly depends whether or not the interval is indeed long enough. It is apparently a difficult task to make such a decision, especially in real time and adaptively by taking all relevant factors into account. However, the prime aim of conducting this study is to deal with the case with short but intermittently occurred satellite deficiency intervals. Only for this kind of case, the performance improvement by employing the proposed SPD would be significant. 2) The complexity of the SPD at satellite deficiency epochs, namely the global RD, is low. In fact it is only a low-dimension recursive least squares problem with 3/2/1 parameters and 6/4/2 measurements. The only difference from an ordinary recursive least-squares estimation is the consideration of the special kind of coloredness of the measurement noises (the time differenced carrier phase measurements only); or to be more specific, only a cross covariance is additionally calculated, stored and update. A final note is that, for an interval without any visible satellites, none of the CSCF techniques can perform any information propagation and historical information has to be abandoned. This is reasonable, because for GNSS-only applications without introducing external aid, on media is available to carry the historical information. With this kind of interval, both conventional RD and the proposed SPD can re-initialize at an epoch when even one satellite becomes visible; however conventional PD has to re-initialize only at an epoch with restored satellite sufficiency. 5. Conclusion CSCF is an effective method to process the two basic types of GNSS signals: code pseudorange and carrier phase. This study aims to improve the positioning performance of GNSS-only CSCF technique under harsh environments with intermittent satellite deficiencies. By GNSSonly we mean that no aid is introduced, either from external sensors or from internal models. By satellite deficiency, we mean that there are visible satellites at some epochs; however the
4946
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
number is less than four. The PD filtering is superior to the RD one in terms of accuracy and insensitivity of satellite changes; however, the former is not applicable at satellite deficiency epochs. So the positioning performance improvement in this study is achieved by employing the generally superior PD filtering, however with special modifications tailored to consider the intermittently occurred satellite deficiencies. This modification is a seemingly simple switching, namely from the PD filtering for sufficiency epoch to the RD filtering for deficiency epoch and conversely. However, in the proposed switching method, the RD filtering for the deficiency epochs is a globally optimal one with all pseudoranges treated together as a vector in contrary to the conventional single-channel RD filtering which is only locally optimal. This global RD filter is necessary to preserve the correlations among different channels. This correlation indeed exist as at the first epoch of a deficiency interval, the smoothed pseudoranges is calculated using the same parameter vector estimate at the previous epoch. Real BDS data is processed and the superiority of the proposed method is validated, namely with higher steady state and transient accuracies than conventional RD and PD filtering, respectively. Acknowledgments This research is supported by the Fundamental Research Funds for the Central Universities of China under grant 2018QNA22. References [1] Z. Zhou, B. Li, Optimal Doppler-aided smoothing strategy for GNSS navigation, GPS Solut. 21 (2017) 197–210. [2] R.R. Hatch, The synergism of GPS code and carrier measurements, in: Third International Geodetic Symposium on Satellite Doppler Positioning, New Mexico, 1982, pp. 1213–1232. [3] P.J.G. Teunissen, The GPS phase-adjusted pseudorange, 2nd international workshop on high precision navigation, Stuttgart, 1991, pp. 115–125. [4] C. Gunther, P. Henkel, Reduced-noise ionosphere-free carrier smoothed code, IEEE Trans. Aerosp. Electron. Syst. 46 (2010) 323–334. [5] G. Chang, T. Xu, Y. Yao, Q. Wang, Adaptive Kalman filter based on variance component estimation for the prediction of ionospheric delay in aiding the cycle slip repair of GNSS triple-frequency signals, J. Geod. 92 (2018) 1241–1253. [6] P.Y.C. Hwang, G.A. McGraw, J.R. Bader, Enhanced differential GPS carrier-smoothed code processing using dual-frequency measurements, Navigation 46 (1999) 127–137. [7] G.A. McGraw, Generalized divergence-free carrier smoothing with applications to dual frequency differential GPS, Navigation 56 (2009) 115–122. [8] B. Park, C. Lim, Y. Yun, E. Kim, C. Kee, Optimal Divergence-free hatch filter for GNSS single-frequency measurement, Sensors 17 (2017) 448. [9] J.Y. Lee, H.S. Kim, K.H. Choi, J. Cho, H.K. Lee, Adaptive position-domain carrier-smoothed code filter based on innovation sequences, IET Radar, Sonar Navig. 8 (2013) 336–343. [10] H.K. Lee, C. Rizos, Position-domain hatch filter for kinematic differential GPS/GNSS, IEEE Trans. Aerosp. Electron. Syst. 44 (2008) 30–40. [11] S.B. Bisnath, R.B. Langley, Precise orbit determination of low earth orbiters with GPS point positioning, in: 2001 National Technical Meeting of the Institute of Navigation, Long Beach, CA, 2001, pp. 725–733. [12] H.K. Lee, C. Rizos, G.-I. Jee, Design of kinematic DGNSS filters with consistent error covariance information, IEE Proc.-Radar, Sonar Navig. 151 (2004) 382–388. [13] H.K. Lee, C. Rizos, G.-I. Jee, Position domain filtering and range domain filtering for carrier-smoothed-code DGNSS: an analytical comparison, IEE Proc.-Radar, Sonar Navig. 152 (2005) 271–276. [14] W. Liu, X. Wang, Z. Deng, Robust centralized and weighted measurement fusion Kalman estimators for multisensor systems with multiplicative and uncertain-covariance linearly correlated white noises, J. Frankl. Inst. 354 (2017) 1992–2031.
G. Chang, T. Xu and C. Chen et al. / Journal of the Franklin Institute 356 (2019) 4928–4947
4947
[15] L. Jiang, L. Yan, Y. Xia, Q. Guo, M. Fu, L. Li, Distributed fusion in wireless sensor networks based on a novel event-triggered strategy, J. Frankl. Inst. (2018) In Press. [16] G. Wang, N. Li, Y. Zhang, Maximum correntropy unscented Kalman and information filters for non-Gaussian measurement noise, J. Frankl. Inst. 354 (2017) 8659–8677. [17] W. Li, Y. Jia, J. Du, Distributed extended Kalman filter with nonlinear consensus estimate, J. Frankl. Inst. 354 (2017) 7983–7995.