Optimization Algorithm for Kalman Filter Exploiting the Numerical Characteristics of SINS/GPS Integrated Navigation Systems.

Aiming at addressing the problem of high computational cost of the traditional Kalman filter in SINS/GPS, a practical optimization algorithm with offline-derivation and parallel processing methods based on the numerical characteristics of the system is presented in this paper. The algorithm exploits the sparseness and/or symmetry of matrices to simplify the computational procedure. Thus plenty of invalid operations can be avoided by offline derivation using a block matrix technique. For enhanced efficiency, a new parallel computational mechanism is established by subdividing and restructuring calculation processes after analyzing the extracted "useful" data. As a result, the algorithm saves about 90% of the CPU processing time and 66% of the memory usage needed in a classical Kalman filter. Meanwhile, the method as a numerical approach needs no precise-loss transformation/approximation of system modules and the accuracy suffers little in comparison with the filter before computational optimization. Furthermore, since no complicated matrix theories are needed, the algorithm can be easily transplanted into other modified filters as a secondary optimization method to achieve further efficiency.

1. Introduction

In many SINS/GPS integrated navigation systems, common system models cannot conform to the use conditions of classical Kalman filter due to the “colored noise”. To still be able to use the optimal estimation method, it is necessary to expand the system state model [1,2]. However, as the state dimensions increase, filtering computation will rapidly become expensive and unstable, and the so-called “dimension disaster” could even break out. In order to address this problem, scholars have reported various optimization algorithms which can lower the computational costs and/or enhance the numerical robustness. The typical optimization algorithms contain square-root information filtering SRIF [3], U-D decomposition filtering [4], singular value decomposition filtering [5] and their improved versions [6,7,8]. These modified filters, besides greatly reducing the time consumption, theoretically guarantee the positive definiteness of covariance and effectively avoid the numerical divergence, thus they are widely applied in the theoretical design of higher-order systems. Nevertheless, they require the support of relatively complex matrix theory in derivation and may pose certain difficulty in engineering. On the other hand, although these algorithms have decreased computational costs, their complexity is still O(n3). With further expansion of the state dimensions, the above decomposition algorithms, as a kind of general method, are also being challenged on their real-time capability. Recently, researchers have paid much attention to non-commonality methods, which are designed for some specific applications or based on specific models. For instance, aiming at integrated navigation applications, a kind of reduced-order Kalman filter (RDKF) [9,10,11,12] was promoted to ease the computational load. The idea of these filters is to reduce the model dimension by theoretical/engineering methods. As the filter dimension n is vital to computation time, the reduction of the state order will produce a direct benefit in terms of real time. However, the decrease of the state order will also bring partial accuracy damage, therefore, these methods may not be suitable for some high accuracy-demanding applications. Another kind of effective algorithm optimizes the float-point operations mainly by taking advantage of the sparse matrices (e.g., the transition matrix Φ in SINS/GPS) during the filtering computation [13,14,15]. In these algorithms, the so-called matrix accumulative method (MAM) or other online methods are used and the algorithm complexity can be simplified to O(s) (s/u represents the numbers of nonzero/1 elements in Φ). As the matrices in the high order model normally contain substantial numbers of zero elements, the computational time can be curtailed to an ideal level by usage of this sparse-matrix-based method. Meanwhile, being different from RDKF, the sparse-matrix-based methods do not change the system model and need little engineering approximation, thus they perform better in the accuracy-control aspect. But these methods also have their own limitations:

Though they avoid massive unnecessary multiplying-zero operations, the methods still need extra O(n3) times estimation of zero elements in each matrix multiplication;

Methods based on MAM, a kind of online method, could not seek any deeper optimization online and their efficiency is totally dependent on the number of zeros in the real-time matrix. To enhance the optimization efficiency, the usual way is to set more zero elements by approximation methods, but this comes at the cost of a certain accuracy loss.

In response to the above problems, we present here a new highly-efficient, accuracy-lossless, and engineering-tractable optimization algorithm based on offline derivation and a parallel method exploiting the numerical characteristics of SINS/GPS. In comparison with other algorithms, the proposed method offers numerous advantages: (i) in comparison with the general algorithms, our proposed method requires little complex deduction and is easily understandable; (ii) in comparison with the RDKF optimization method, the proposed method needs little engineering approximation and would be more accurate; and (iii) in comparison with the filter based on MAM, it is free of zero-estimation operations and can deduce some stronger conclusions, thus it is more efficient. Emphasizing on the engineering simplicity, we will introduce our optimization method revolving around the classical Kalman filter and a loosely-coupled model of SINS/GPS. It is apprised here that the chosen model does not imply any constraints on the application scope of our method. Actually, the offline and parallel method, as a pure numerical optimization method, is equally applicable to the extended filter and other complex SINS/GPS models if some similar research on the special model is done. On the other hand, since there is no transformation on any models or equations in our method (the two distinct differences with the normal filtering computations are: (i) unnecessary floating-point operations are directly avoided; and (ii) de-coupled operations are processed in parallel), the precision and robustness of our method would be equal to the classical Kalman filter in theory. Owing to this fact, the efficiency rather than robustness or accuracy is focused on in our discussion later. Section 2 gives a brief account of closed-loop Kalman filtering and a widely-used loosely-coupled model in SINS/GPS. In Section 3, the details of offline derivation and the parallel method based on a block-matrix technique are described. Section 4 and Section 5 deal with the statistics on computational costs and the filter performance on error estimation.

2. Loosely-Coupled SINS/GPS Model Using the Classical Kalman Filter

Depending on the different ways of integration, SINS/GPS are classified into two modes: (i) loose mode; and (ii) tight mode. Loose mode using GPS output to adjust SINS errors is opted for in this research for its distinctive features of high redundancy and easy realization. The emphasis is placed on a common kind of loose model i.e., position-velocity integration with 18 states and six observations.

2.1. State Equation

By using the indirect method, which treats the estimation variables as parametric errors instead of parameters themselves, the state equation is established in Equation (1): where, represents the state vector, including attitude errors , velocity errors , position errors , and the “colored noise” errors: gyro constant drift , Markov process drift , and accelerometer drift ; denotes the collected system noise vector, including gyro drift white noise , Markov driving white noise , and accelerometer drift white noise ; is state transition matrix and is noise coefficient matrix. Combining references [16,17,18,19,20,21,22] with some modification, the variables of the state equation can be described by Equations (2)–(5):

In matrices (4) and (5), all blocks are 3-order square matrices. Blocks , , and stand for the zero matrix, the identity matrix, and the attitude matrix, respectively. The other blocks are described in detail in [16,17,18,19,20,21,22]. In particular, we highlight here two kinds of nonzero blocks and of which the special single-nonzero-vector and diagonal structure may contribute to deeper optimization:

In addition, the variance of noise vector is described as: where, block is also diagonal.

2.2. Measurement Equation

We adopt a position-velocity integration model and define the filter observation vector [16,22]:

Then, we define the measurement equation as follows: where represents the measurement noise vector; denotes the corresponding variance matrix, and is the coefficient matrix. , and are described as follows according to [16,17,18,19,20,21,22]: where:

2.3. Closed-Loop Kalman Filtering

In order to fulfill the requirement of closed-loop Kalman filtering, the state equation (Equation (1)) and measurement equation (Equation (9)) are changed to the discrete forms: where, the transition matrix is converted from ; the noise driving matrix is converted from ; the coefficient is converted from ; and the control vector is added for real-time calibration on the parameters of SINS.

After discretization, the update of is done recursively as follows: where, “~” represents the parameter needed estimation; represents the value before estimation; represents the value after estimation; represents the value after calibration with the vector ; represents real observation vector; and represents Kalman filtering gain.

To minimize the filtering error, state vector is permanently set to zero:

Substitution of Equation (17) into Equation (16), yields:

Apparently, becomes the final output of the Kalman filter with control vector. As is almost directly obtained by measurement, appears to be the key of and even the filter. According to the Kalman filtering theory, can be recursively computed as follows: where is the covariance; and are the discretized forms of and respectively. In fact, Equations (18)–(21) constitute the main process of the closed-loop Kalman filter. This process is generally subdivided into two processes: time propagation Equation (19) and measurement updating Equations (18), (20) and (21). In summary, the process flow of Kalman filter in SINS/GPS is depicted in Figure 1. Since most matrices involved in Figure 1 or Equations (19)–(21) are either 18 × 18 or 18 × 6 matrices, there will be a large-scale and time-consuming computation in each recursive cycle. Therefore, a certain optimization on the computation would be essential to the real-time properties of SINS/GPS.

Figure 1

Iterative closed-loop Kalman filtering process.

3. Computational Optimization

Consistent with Figure 1, optimization are sequentially implemented in the updating processes of , , , , , where the most computational time are consumed.

3.1. Compute

According to [16,17,18,19,20,21,22], we can calculate as follows:

Here, the 3rd order Taloy equation is selected for mainly two reasons:

(1) Equation (22) meets the requirements of most mid-high-accuracy systems;

(2) Equation (22) can deduce the universal optimization conclusion.

It’s not difficult to prove that would have a uniform numerical structure independent from which of the 3rd and higher order equation is chosen. On the other hand, lower order equations are obviously a special case of 3rd order equation. Therefore, the conclusion based on 3rd equation is suitable for the model based on the different order of equations too.

As the computation of Equation (22) is centered on the 18 × 18 matrix , it is necessary to inspect the numerical characteristics of before computational optimization: matrix , in which only 64 of the total 324 elements are non-zero, is clearly a sparse matrix. A substantial reduction of useless multiplying-zero operations can be achieved by directly expanding and operating every element in deduction (namely “direct derivation”), but because of the high order, direct derivation is quite complicated to handle and is not advisable in practice. Fortunately, the 3-order-block form of , having 23 zero blocks among the total of 36 blocks (see Equation (4)), is still a sparse matrix structurally. Therefore, we can carry out our offline optimization using an indirect method, in which the blocks instead of the elements are treated as the atomic unit in deduction. In this way, we can largely decrease the derivation complexity and retain most of the efficiency of a direct method at the same time. For these reasons, the indirect method is applied to all the derivation processes in this paper. In addition, the nonzero blocks and (see Equation (6)), which need few operations in multiplication by any 3-order block (no more than 32 float-pointing multiplications and 2 × 3 additions), are also introduced to the derivation for higher efficiency.

It is a fact that manual derivation with indirect method is still a cumbersome job in view of multiple matrix multiplications (e.g., , ), so we handle our derivation by a computer program instead. With the powerful symbolic-computation function of some well-known math software, e.g., MATLAB and Maple, the complicated derivation can be easily accomplished. After necessary programming, the final derivation results are shown in Equations (23) and (24):

Also, the computation of blocks and is described in a computerized sequence as follows: whereas all products involving block , or are zero blocks and need no real-time computations; products involving the special blocks, e.g., , ( single-nonzero –vector block), , , or (diagonal blocks), are simplified and are partly indicated in symbol .

Substituting of Equation (4), Equation (23), and Equation (24) in Equation (22), yields: whereas:

Here, and are single-nonzero-column blocks; is the sum of and single-nonzero-column block; and are diagonal blocks. These blocks would be specially used in other processes.

In conclusion, with offline derivation based on 3rd order block method, nearly half amount of blocks is proved to be constant blocks: or (see Equation (26)), thus, the computations on them are simply eliminated. Additionally, by using the properties of special blocks , , , , , and introducing the intermediate variables , , real-time computations on the other half can be further reduced.

3.2. Compute

According to [16,17,18,19,20,21,22], the 2nd order approximate of are computed as follows: where:

Obviously, and are both symmetric. Thus, for matrix , only its upper (or lower) triangle elements need to be computed in real time. Substitution of (Equation (29)) in (Equation (28)), yields: where:

In Equation (31), and are apparently diagonal blocks that may be used to simplify other processes. As stated above, we optimize the computation of mainly by applying some special numerical properties to the offline derivation (e.g., the sparsity of , the simplicity of diagonal blocks , , , , and the symmetry of ) and much of the computation time is saved.

3.3. Compute

For simplicity, and will be written in a unified form (this is allowed because of that and are actually the same covariance matrices in different state and share the same memories). Similar to or , is written in 3rd-order-block-matrix form:

According to its definition, is a symmetric matrix:

Thus, only the lower (or upper) triangle part needs substantive computation. Regardless of at first, Equation (19) is written as:

Substitution of Equation (26) and (32), yields: where the symbol represents intermediate items that would be temporarily stored in memory to avoid repeated computation; the symbol represents the repeated item that have been computed before and are available by directly accessing and the symbol indicates that the special properties of are exploited to perform further optimization.

After Equation (35) is done, should be added. Substitution of Equations (30) and (32) in Equation (19), yields:

Above all, the computation of is divided into two steps: Equations (34) and (36). As the latter has few computations, optimization mainly focuses on the former. Employing the sparsity of , the symmetry of and the repeatability of intermediate items, the real time of calculating Equation (34) is greatly improved.

3.4. Compute

Substitution of Equations (11), (12) and (32) in the inversion part of Equation (20), yields:

Since has a low order, few computations are needed during inversion (the complexity is that can be neglected in comparison with other procedures) and the details of inversion are not discussed accordingly. Here, we use the symbol to represent the result of inversion directly:

Matrix is easily proved to be symmetric and can be written in the 3-order-block-matrix form as:

Then can be computed as:

Define as: then:

The derivation result shows that, by using the block-matrix optimization technique, only a simple process (Equation (37), (38) and (42)) is needed in updating. Compared with the complicated Equation (20) which needs multiple 18-order matrix production, our method seems much more simple and efficient in computing.

3.5. Compute

Transform Equation (21) into:

Obviously, the pressure in computation mainly comes from . Fortunately, is a symmetrical matrix (this is an evident inference according to Equation (43) in the condition that , are both symmetrical) so that only its upper (or lower) triangular part needs real-time updating. Define as the variance’s increment :

Substitution of Equations (12), (32) and (41), yields: then:

In conclusion, taking advantage of the symmetry of and the sparsity of , we deduct a simple form of computing : Equations (45) and (46). Equation (46) is apparently costless as only the matrix subtraction is involved. Equation (45) involves matrix products, but requires no computations on most blocks of (except the blocks in 2nd and 3rd rows) and proves to be highly efficient.

3.6. Parallel Computation

The optimization results (see Table 2, Table 3, Table 4, Table 5 and Table 6) reveal that, computation on costs nearly half of the processing time and becomes a bottleneck for further efficiency. The low efficiency mainly arises from that becomes increasingly dense and can provide few zero blocks after initial filtering cycles. In low-accuracy applications, some researchers may force the some matrix elements to be zero for low computational costs with the engineering approximate method. However, to some accuracy-demanded applications, methods with more accuracy and efficiency are still looked for. Parallel computation is one of the feasible approaches.

Table 2

Optimization efficiency of updating process.

	Offline-Derivation Method	Traditional Method
Multiplication	42 × 3^3 + 36 × 3^2 + 13 × 3 × 3^2=1809 (I), 10%	3 × 18^3 + 3 × 18^2 = 18468 (II)
Addition	42 × 2 × 3^2+9 × 2 × 3+13 × 2 × 3^2=1044 (I), 6%	3 × 17 × 18^2 + 3 × 18^2 = 17496 (II)
Memory	32 × 3^2 × 8 Byte=2.3 KB (III)	3 × 18^2 × 8 Byte =7.6 KB (III)

Table 3

Optimization efficiency of updating process.

	Offline-Derivation Method	Traditional Method
Multiplication	24×3^3+18×3^2+12×3^2=918, 6%	2×18^3+18^2+18×6^2+6×18^2=14580
Addition	24×2×3^2+9×2×3+13×3^2=603, 4%	2×17×18^2+18^2+5×6×18+17×18×6=13716
Memory	13×3^2×8 Byte=0.9 KB	3 × 18^2 × 8 Byte=7.6 KB

Table 4

Optimization efficiency of updating process.

	Offline-Derivation Method	Traditional Method
Multiplication	120×3^3+37×3^2=3573, 31%	2×18^3=11664
Addition	120×2×3^2+18×2×3+134×3^2=3474, 31%	2×17×18^2+18^2=11340
Memory	60×3^2×8 Byte=4.2 KB	3×18^2×8 Byte=7.6 KB

Table 5

Optimization efficiency of updating process.

	Offline-Derivation Method	Traditional Method
Multiplication	24 × 3^3 + 6 × 3^2 = 702, 11%	18 × 18 × 6 + 18 × 6 × 6 + 6 × 18 × 18 + 6 × 18 × 6 = 5184
Addition	24 × 2 × 3^2 + 12 × 3^2 = 540, 11%	17 × 18 × 6 + 17 × 6 × 6 + 5 × 18 × 6 + 17 × 18 × 6=4824
Memory	16 × 3^2 × 8 Byte= 1.1 KB	（3 × 6 × 18+6 × 6）× 8 Byte = 2.8 KB

Table 6

Optimization efficiency of updating process.

	Offline-Derivation Method	Traditional Method
Multiplication	42×3^3+6×3^2=1188, 15%	6×18×18+18^3=7776
Addition	42×2×3^2+42×3^2=1134, 15%	5×18^2+17×18^2+18^2=7452
Memory	16×3^2×8 Byte=1.1 KB	（3×6×18+6×6）×8 Byte=2.8 KB

Before discussing the parallel method, we perform an inspection of the whole computation processes in order to clarify the dependency of every parameter on . The numerical relationship between and other parameters can be consulted in Table 1.

Table 1

Dependency of every real-time parameter on blocks.

ϕ_ n, Q_n, P_n(−)	K_n (Equations (40) and (41))	P_n(+_e) (Equation (43))	ΔP_ij (Equation (45))
Irrelevant to P_n,ij	Dependent on P_n,i2, P_n,i3 only	Dependent on P_n,ij	Dependent on P_n,2i, P_n,3i only

Except for , all parameters are either irrelevant to or only relevant to 2nd and 3rd block rows (columns) of . Taking this fact into account, we can call 2nd and 3rd block rows (columns) “useful” data and treat the rest as “useless” (it is emphasized here that “useless” only refers to the requirement of most computation processes, but not that the data is really useless and needs no updates).

Furthermore, because is a symmetrical matrix, the “useful” blocks are only 2nd, 3rd block columns in fact. If “ ” represents “useless” blocks, can be written as:

In measurement updating (Equations (18), (20) and (21)), all parameters except are decoupled with the “useless” data (Equation (19)) which belong to the time propagation. Moreover, to the coupled its time-consuming part is also independent of the “useless” blocks. In other words, with some appropriate modification of measurement updating and time propagation, we can produce a new decoupling mechanism of the two procedures. The modification contains three steps:

Step1. Subdivide time propagation.

Classify the blocks of into two kinds: “useful” data and “useless” data;

Step 2. Subdivide measurement updating.

Subdivide into two processes: computing (Equation (45)) and adding (Equation (46));

Step 3. Restructure measurement updating.

Restructure measurement updating process as (instead of ) (Equations (18) and (21)).

With the above steps, the restructured process and the “useless” data updating process are completely numerically decoupled and can be computed in parallel. In this way, the efficiency bottleneck in updating is broken.

Furthermore, in order to achieve the parallel computation on the “useful” data as well, we similarly subdivide the process of updating “useful” data into two steps: computing (Equation (34)) and adding (Equation (36)). Apparently, computing and updating (Equation (31), which belongs to time propagation) can be processed in parallel. In summary, the process of Kalman filtering with parallel computation is described as Figure 2.

Figure 2

Filtering process with parallel computation.

It is pointed out here that the parallel method in this paper has a fundamental difference with some other parallel methods [23,24]. The others owe their decoupling to the so-called “mandatory delay”, which may cause certain accuracy damage. On the contrary, the decoupling in this paper needs neither “mandatory delay” nor engineering approximation and is achieved mainly by introducing the “useful” data and subdividing the computation process. Essentially, it is the numerical characteristic of SINS/GPS that brings up this special decoupling. Therefore, the parallel processing in this paper is an accuracy-lossless method.

In respect of efficiency, the parallel optimization method can reduce much execution time needed in updating both and the “useless” blocks (including 918 multiplications and 603 additions in updating , and 1089 multiplications and 1107 additions in updating “useless” blocks). In comparison with non-parallel methods, efficiency of parallel method increases by about 25%.

4. Optimization Efficiency

According to the derivation result, we can figure out the optimization efficiency by counting and comparing the computational counts needed in both the proposed method and the existing algorithm. In detail, the discussion about counts proceeds on each parameter updating procedure, respectively. As matrix multiplication is mainly involved, the calculation will focus on the floating-point operations of multiplication and addition. Table 2, Table 3, Table 4, Table 5 and Table 6 will show the statistical details of each parameter. Then, Table 7 and Table 8 give the global efficiency and the comparison with other general algorithms, respectively. As shown in the tables, the offline-derivation and parallel method saves about 90% of computational time and 66% of memory space, while keeping the required accuracy level. Compared with general decomposition optimization algorithms, the proposed method needs no complicated matrix theory but producing higher efficiency. Under the condition that high hardware precision has guaranteed the numerical robustness, our algorithm appears quite simple and practical.

Table 7

Optimization efficiency of offline-derivation and parallel method.

	Offline-Derivation & Parallel Method	Traditional Method
Multiplication	1809+2484+702+1188=6183 (IV), 10.7%	18468+14580+11664+5184+7776=57672
Addition	1044+2475+540+1134=5193 (IV), 9.5%	17496+13716+11340+4824+7452=54828
Memory	2.3+0.9+4.2+1.1+1.1=9.6 KB	7.6+7.6+7.6+2.8+2.8=28.4 KB

In (IV), 2484 (2475) is the operation times of “useful” blocks updating.

Table 8

Comparison of offline-derivation & parallel method and general optimization filters.

	Multiplication	Addition
SRIF filtering	76×18^3+36×18^2=18468 (V)	76×18^3+36×18^2=18468 (V)
U-D decomposing filtering	12×18^3+12×18^2+10×18^2+8×18=6462 (VI)	12×18^3+12×18^2+9×18^2+9×18=6156 (VI)
SVD filtering	26×18^3+78×18^2=176904 (VII)	26×18^3+78×18^2=176904 (VII)
offline-derivation and parallel method	2484+702+1188=4374 (VIII)	2475+540+1134=4149 (VIII)

In (I), 42 × 33 (42 × 2 × 32), 36 × 32 (9 × 2 × 3) and 13 × 3 × 32 (13 × 2 × 32) represent the counts on multiplication (addition) of items in Equation (25) without special blocks, of items in Equation (25) with special block and of items in Equation (27) respectively, and only 32 & 33 items are counted; In (II), the first item represents the counts on operations of 18-order-matrix multiplication, the second item represents that of multiplication by scalars or matrix addition; In (III), data are assumed to be double-precision (8-byte word length); the percentages represent the efficiency in comparison with the traditional method.

(V) is explained in Reference [3], while (VI) in Reference [25], Reference [26] and (VII) in Reference [5]; Point out here that, (VI) relies on the precondition that is already an upper-triangle matrix in 9-order-block form (what is the conclusion of offline derivation), otherwise, 9702 multiplications and 9396 additions are actually needed; concerning (VIII), the procedure of computing is no longer counted as no optimizations are against in the general algorithms.

5. Simulation

To evaluate the method validity, we design a practical loosely-coupled SINS/GPS program in CCS v5.5 (a well-known programming tool on DSP platform). The program comprises three modules:

Sensor data sampling module.

SINS algorithm module. This module is designed to implement the calculation of attitude, velocity and position independently (only the inertial sensors data is needed).

Kalman filter module. This part is for data fusion. The system module is designed according to Section 2, and the computation process is programmed using the offline and parallel method described in Section 3. Besides, for accuracy evaluation, the classical KF is also designed as a control group.

With the developed program, we can simulate filtering process and evaluate the performance of our method in CCS. The simulation parameters and conditions are set as follows:

DSP TMS320C6722, a famous float-point CPU, is chosen in CCS as the computing device.

The carrier of the navigation system is assumed to be static. With this assumption, the gyro/accelerometer would detect no valid rotate/velocity rate aside from devices noise. Therefore, noise is considered the sampling data driving the SINS/GPS program. As the carrier is static, velocity/position should stay at the zero/initial value in theory, thus, this assumption can largely facilitate the evaluation of program results.

Velocity/position measured by GPS is constantly set: Velocity/position calculated by SINS is initialized:

IMU data rate (noise frequency): 100 Hz; SINS velocity/position updating rate: 50/20 Hz; filtering rate: 10 Hz.

With the above setting, we run the program in CCS for 100 s and record the navigation outputs (e.g., the X-axis velocity and latitude). The simulation results are illustrated in Figure 3 and Figure 4.

Figure 3

Latitude calculation by different methods.

Figure 4

Velocity calculation by different methods.

The results show that the proposed method effectively restrains the unbounded error accumulation in SINS and keep the estimated variables around their mean values. On the aspect of estimation performance, our method is nearly at par to the classical KF (the proposed-method curve does not completely overlap the classical-KF one because of less calculation and less truncation errors). This conclusion agrees with the theoretical analysis: our numerical method, needing no modification of the system module and no engineering approximation, causes little damage to the estimation accuracy, but from the point of view of real-time computation, our method is much more efficient than the classical KF. In the DSP program, the proposed method needs 64,820 system clocks in each filtering cycle while classical KF needs 390,480 clocks. If a clock frequency of 200 MHz, the typical frequency of a DSP TMS320C6722 is chosen, then our method will cost only 324 μs in each filtering operation. Such a level of real-time efficiency is highly valuable to the many accurate navigation applications.

6. Conclusions

Starting from a common classical-Kalman-fitler-based SINS/GPS model (position-speed integration with 18 states and six measurements), we present an optimization algorithm based on the system’s numerical characteristics. The algorithm employs a block-matrix technique in offline derivation, where special blocks (e.g., zero blocks, diagonal blocks, etc.) are used to simplify the calculation. In this way, plenty of invalid multiplying by zero and repeated operations are avoided offline. Furthermore, aiming at the time-consuming update of , a novel parallel method is implemented by defining “useful” data and subdividing computational process, and the entire filtering procedure is greatly accelerated. The statistical analysis and simulation results reveal that the offline algorithm, coupled with the parallel method. can reduce the CPU processing time by 90% and the memory usage by 66%. Compared with several general decomposition-optimization algorithms, the proposed method requires no complex matrix theory rather more efficiency. The distinguished feature of the algorithm is that no modification of the system module and no engineering approximation are needed, thus it causes little harm to the base filter. It is pointed out that although the derivation is based on a specific integration model, the algorithm as a complete numerical optimization approach can be transplanted to other advanced models of SINS/GPS or other extended filters. With the powerful symbolic-operation function of the MATLAB program, researchers can manipulate formulas involving high order matrices in an easy way. Consequently, the proposed method is an engineering-tractable approach with high efficiency and high precision for SINS/GPS integrated navigation systems.