Reducing Computational Load in Solution Separation for Kalman Filters and an Application to PPP Integrity

Abstract

This paper investigates two techniques to reduce the computational load of running multiple fault tolerant Kalman filters in order to provide integrity. These approaches are then exploited in the implementation of a solution separation integrity monitoring algorithm in a PPP Kalman filter solution. We evaluate the techniques using GNSS data collected in static and driving conditions. In our scenarios, these techniques lead to computational load reductions of at least 70% at the expense of protection level degradations of about 50%. INTRODUCTION Until recently, Precise Point Positioning (PPP) techniques [1] have mostly been used to provide high accuracy. There is a growing interest in translating the benefits of PPP to integrity and enabling its application to safety critical applications in rail, automotive, maritime, and even air navigation [2], [3], [4], [5]. In [5], we demonstrated how techniques developed for aviation applied to PPP can produce meter-level protection levels in automotive and aviation scenarios. This was achieved by implementing an integrity monitoring algorithm based on solution separation, akin to the one used to analyze Advanced RAIM performance [6], to the PPP Kalman filter solution. The principle of solution separation is to run a bank of filters, where each filter is fault tolerant to a fault or set of faults. The fault detection statistic is the difference between each of these solutions and the all-in-view solution. In addition to their optimality properties [7], solution separation algorithms offer a straightforward proof of integrity, and good performance [5]. However, they can also be expensive in terms of memory and processing time, because they require the receiver to compute a bank of filters (or a process computationally equivalent to a bank of filters, as in [11]). In the worst case, the computational load will be proportional to the number of filters. In [5] we showed that it was possible to dramatically reduce the cost of running the bank of filters: depending on the filter complexity (that is, the number of estimated states), we could run 20 to 50 additional filters for the cost of one. This was obtained by exploiting the fact that, in PPP, many of the elements in the computation (error models, corrections, etc) are common to the all filters, so that it is sufficient to compute them once for the all-in-view filter. Also, all the measurements are linearized with respect to the all-in-view position solution, which further simplifies the subset solution filters. The goal of this paper is to introduce and investigate techniques to reduce even further the cost of the solution separation for Kalman filter solutions. When the number of states is large (larger than 50), which is the expectation in a PPP multifrequency user algorithm, there are at least two steps that are computationally expensive: the determination of the Kalman gain, and the determination of the new error covariance. The first technique under investigation consists in using a set of suboptimal filters for the subset solutions (instead of the optimal filter), where the Kalman gain for each subfilter is derived from the all-in-view and does not require a full matrix inversion. The second technique that we will evaluate is the consolidation of faults into a few subsets. This is another type of suboptimal subset solution that has already been evaluated for snapshot solutions in Advanced RAIM [8], [9]. For the experimental evaluation, we will use our sequential PPP filter implementation [5], which is based on a simple extended Kalman filter with estimated parameters comprising the receiver position, clock biases for each constellation in use, a tropospheric delay, and float ambiguities for each tracked carrier phase. Dual-frequency measurements are incorporated from GPS, GLONASS, Galileo, and BeiDou. The precise orbit and clock estimates are drawn from the IGS MGEX analysis centers. SUBOPTIMAL SUBSET SOLUTIONS: FIRST APPROACH For the approach outlined here, we only consider the measurement update step of the Kalman filter. For the all-in-view filter (indicated by the index 0), we have the following Kalman filter equations:           0 0 0 0 1| 1 1| 1| 1 1 1| ˆ ˆ ˆ T t t t t t t t t t x x C G W y Gx           (1) Where   0 1| 1 ˆ t t x   is the a posteriori state estimate,   0 1| ˆ t t x  is the a priori estimate, G is the observation matrix, W is the inverse of the measurement noise matrix, y is the vector of measurements, and   0 1| 1 t t C   is the error covariance of the a posteriori state estimate. We have:         1 1 0 0 1| 1 1| T t t t t C C G WG        (2) Where     1 0 1| t t C   is the error covariance of the a priori estimate. Subset filter solutions In this paper we will consider probability faults of 10 and 10 per hour. For an integrity of 10 per hour, this means that we need to compute the solution separation statistic of all one out subsets in the case of 10 and two out subsets for 10 [8]. The subset filter Kalman filter equations (indexed by k) are similar to the all-in-view ones:                 1| 1 1| 1| 1 1 1| ˆ ˆ ˆ k k k k T k k k t t t t t t t t t x x C G W y G x           (3)               1 1 1| 1 1| k k k k k T t t t t C C G W G        (4) The only difference with respect to the all-in-view is that we only use a subset of the available measurements to update the state estimate. One of the most onerous steps in this process is the computation of the covariance as written in Equation (4) . As can be seen, we need at least two matrix inversions, where both matrices are n by n with n being close to a 100. (We note that the geometry matrix is often as large, so a matrix update formula will not substantially reduce the computational load). The first approach consists in using a suboptimal filter for   1| 1 ˆ k t t x   instead of the optimal one defined above. More precisely, we define it as follows:                   1| 1 1| 1| 1 1| ˆ ˆ ˆ k k k k T k k k k t t t t t t t t x x G W y G x          (5) Where the matrix   1| 1 k t t    is no longer given by Equation (4). Instead, we attempt to find a matrix that will result in a reasonable estimator but that is cheaper to compute. One possible approach is to compute this matrix as if the prior of the estimated states was given by the prior of the all-in-view filter:               1 1 0 1| 1 1| k k k k T t t t t C G W G         (6) The advantage is that this matrix can be obtained without a full matrix inversion. We have:               1 1 0 1| 1 1| k k k k T T T t t t t C G WG G W G G WG           (7) In most cases, the rank of the matrix       k k k T T G W G G WG  is considerably smaller than the rank of     1 0 1| T t t C G WG    . For example, in the one out case of our PPP filter, the rank of this matrix is 4. We can write:             k k k k k k T T T G WG G W G G W G   (8)               1 1 0 1| 1 1| k k k k T T t t t t C G WG G W G          (9) Using the Woodbury matrix identity, we get:                       1 0 0 1 0 0 1| 1 1|t 1 1|t 1 1|t 1 1|t 1 k k T k k k k T t t t t t t C C G W G C G G C                 (10) The use of this formula can speed up the calculation of   1| 1 k t t    , because the matrix to invert is usually much smaller than the whole covariance matrix. Standard matrix inversion algorithms require around 2/3n basic operations, so the computational load is significantly reduced (from almost a million to less than a hundred). The new Kalman gain is given by:                                     0 0 0 1|t 1 1|t 1 1| 1 0 1| 1 1| 1 1 t t k T k k T k k T k t k k T k k k k t k T t t t t C G W C G K G W I W G C G G W C G W                 (11) Where we have highlighted what has already been computed in the all-in-view filter. New subset covariance As opposed to the optimal filter, the covariance after the update is not given by   1| 1 k t t    . Instead, it is given by:           1 1| 1 1| T k k k T t t t t C I K C I K KW K         (12) SUBOPTIMAL SUBSET SOLUTIONS: SECOND APPROACH The second approach can be described much more succinctly. It consists in grouping the faults so that we do not need to run as many filters. For example, instead of running a filter for a fault in satellite i and another one for a fault in satellite j, we run a filter that is fault tolerant to both i and j. This will result in a weaker solution position, and therefore larger protection levels. This second approach can be considered a suboptimal subset solution approach because each fault is accounted by a suboptimal filter. For this paper, the groups were formed based on the PRN number, which is mostly equivalent to a random grouping with regard to geometry. DATA AND PROTECTION LEVEL CALCULATION We used two types of GNSS data: one collected by a static receiver and one collected by a receiver installed in a car. The GNSS data collected in road conditions is described in [5] and briefly summarized here: • Receiver: NovAtel OEM 7500 • 1 Hour Driving Data on March 1, 2018 • GPS (L1 C/A -L2P semi-codeless), GLONASS (L1 C/A-L2P) at 1 Hz • Truth positions provided by NovAtel OEM729 with tactical-grade IMU with forward and reverse processing Specifically, we choose the open sky conditions. The static GNSS corresponded to the following conditions: • Receiver: Trimble NetR9 • 6 hours of static data on November 7, 2018 at Stanford • GPS (L1C-L2W), GLONASS (L1C-L2P) at 1 Hz • Truth position from IGS station solutions The error models are and the protection level calculation is also described in [5]. It is a straightforward adaptation of the ARAIM algorithm described in [6] to a Kalman filter solution. The algorithms were r