Reduced Navigation Error Using a Multi-Sensor Fusion Technique and Its Application in Synthetic Aperture Radar

Modern navigation solutions rely on a combination of inertial measurement units (IMUs) and a global navigation satellite system (GNSS) receiver to estimate a navigating body's position. In order to produce a high-fidelity solution, the current approach is to utilize a single ultra-low bias navigational grade IMU in the system. However, these high-quality IMUs are expensive, bulky, heavy, and require significant power consumption. This article proposes the fusion of multiple lower-quality IMUs to achieve near-identical or better positional accuracy as a single high-quality sensor to minimize cost, size, weight, and power (C-SWaP) without sacrificing the positional estimation accuracy. The primary focus is on a generalized method to fuse multiple position estimations from multiple co-located IMUs for a single navigating body. The proposed fusion algorithm is applied to simulated data produced by three precision micro-electromechanical systems (MEMS) grade IMU modules (Analog Devices ADIS16465) from two different simulated flight paths. The absolute error is calculated between the position estimation generated by the proposed algorithm and the “truth” position provided by the simulation to determine the accuracy of the final result. The error of the proposed algorithm using the Analog Devices modules is then compared to the error between a single navigational grade IMU (NovAtel IMU-ISA-100 C) and the simulated “truth” position. The results show that using only three precision MEMS grade IMUs; the proposed method can produce identically accurate position estimations as a navigational grade IMU while drastically reducing C-SWaP. This result is further validated in measured data from an instrumented test setup using the above mentioned IMU configurations fused with a standard GPS and compared to a real-time kinematic GNSS setup used as a the third-party ground truth. In addition, the proposed method is further validated by integrating the two navigation systems with a Ku-band radar to produce synthetic aperture radar images. The image produced using the multi-IMU configuration as opposed to the navigation-grade IMU is more focused. Given these results, the proposed method has been validated in theory, simulation, and measurement, resulting in an order of magnitude reduction in cost, size, and power consumption, as well as a three-and-a-half times weight reduction of the overall IMU solution.


I. INTRODUCTION
Modern navigation popularly uses the global navigation satellite system (GNSS) to determine timing and the location of a navigating body.However, GNSS traditionally updates at 1 Hz, oftentimes referred to as one pulse-per-second (PPS), and even more advanced real-time kinematic (RTK) GNSS systems are limited to low update frequency rates on the order of 5-50 Hz.Moreover, GNSS requires an open sky condition to obtain a reliable signal.Both of these conditions are the primary limitation of GNSS [1], [2].To overcome these limitations, an inertial measurement unit (IMU), commonly used in conjunction with GNSS, provides data at a higher frequency, typically on the order of 100 to 200 Hz.The IMU data is used to measure acceleration and orientation which is used with the latest GNSS solution to estimate the navigating body's current position (latitude, longitude, altitude), and the faster update rate of the IMU data allows for position estimates to be reported more often than GNSS.However, IMUs will drift or "walk," and the associated errors are compounded as the position estimate is made over time [3], [4], [5].Moreover, the position is calculated by integrating the acceleration measurements provided by the IMU, and as a result, the positional error will increase at a significant rate [6].Therefore, data from the GNSS and IMU are used in conjunction or fused to obtain a more precise position estimate [7], which has been shown to be a viable solution in autonomous navigation [8], [9].
This data fusion concept is commonly used to combine the data from the GNSS and an IMU to produce a single-position solution.In this particular scenario, data fusion combines the accuracy of the GNSS solution as well as the faster, more stable IMU measurement resulting in a single position solution that is stable, accurate, and updated rapidly.Traditional IMU and GNSS fusion applications combine the measurements from a single IMU with the GNSS solution.In order to increase the accuracy of the position solution, a more accurate IMU, such as an ultra-low bias navigational grade IMU, can offer superior sensor performance.While a navigational grade IMU will produce the best positional estimation, these modules are typically large in size, expensive, heavy, and consume significant power.Therefore, there is a direct trade-off between positional accuracy and an increase in cost, size, weight, and power (C-SWaP).This trade-off becomes a significant challenge for navigating bodies, such as unmanned aerial vehicles, that may require highly accurate position solutions but are SWaP-constrained.
This article proposes a generalized multi-IMU particle filter fusion algorithm based on sequential importance resampling (SIR) for fine-resolution position estimation.This position estimation solution relies on the data fusion of multiple IMUs rather than the positional estimation of a single IMU in the classic "brute force" method.Thus, lower quality micro-electromechanical system (MEMS) grade IMUs with increased bias errors can be used, which drastically reduces C-SWaP, but still allows for a highly accurate position estimation to be made after fusing the data to extract the effective truth location.The particle filter is chosen as the preferred fusion mechanism as it portrays excellent performance in a non-linear and non-Gaussian environment.Another significant advantage of the proposed algorithm is that it can work using similar or dissimilar IMUs.
Airborne applications provide complex motion profiles that can provide challenges in the formation of SAR images.In order to capture the complex scenario, simulated flight data is used to verify the proposed algorithm's mathematical implementation and positional accuracy.IMU data for the flight paths are simulated for three co-located Analog Devices precision MEMS grade IMUs.The data is fused with the implemented algorithm, and the fused positional accuracy is calculated compared to the simulated truth values.Furthermore, the error is calculated between the simulated truth position and the position estimation data produced by a navigational grade IMU.
The motivation for the higher positional accuracy is to obtain highly focused images in synthetic aperture radar (SAR) systems.In SAR, the effects of positional errors are evident in the final image due to misalignment of phase from illuminated scatterers, as explained in [10].Traditionally, large GNSS/IMU systems are utilized to obtain a higher positional accuracy.However, in C-SWaP constrained systems, the larger GNSS/IMU systems become problematic.In efforts to miniaturize SAR systems, the technique in this article aims to reduce the C-SWaP of the navigation system for SAR without sacrificing the accuracy required to obtain a highly focused image.
In addition, the NovAtel and the multi-IMU systems can be instrumented in tandem and post-processed with the same radar data to provide a direct comparison between the two navigation solutions.The resulting SAR images obtained in this work are directly compared and directly validate that the multi-IMU solution is capable of providing better accuracy than the single IMU solution while reducing C-SWaP.
The article is structured as follows: Related works are described in Section II, and the position requirements of SAR are explained in Section III.The proposed multi-IMU fusion algorithm and sensor weight calculations are described in Section IV.Section V, provides a method of estimating the number of sensors required to achieve a theoretical desired accuracy.Section VI introduces the sensor model used for the proposed algorithm.Section VII provides a method of estimating the number of sensors required to achieve a theoretical desired accuracy without assuming the sensor errors are Gaussian distributed.The results and analysis of the proposed method as exercised in simulations are described in Section VIII.Section IX presents the result and analysis of the instrumented setup and the experiment used to compare the proposed technique to the traditional navigation solution.Section X summarizes the methods used to calculate the theoretical performance of the proposed fusion method and verification of the theoretical results through simulated and instrumented results.It should be stated that Sections V through IX reuses some content from the authors dissertation [11] with permission.

II. RELATED WORKS
Complex sensor systems utilize disparate sources of information to produce a more accurate or stable measurement.These systems often integrate measurements from multiple sensors to more accurately measure environmental properties such as location, target tracking, or obstacle detection.For example, the application in [12] utilizes a technique of integrating multiple sensors to produce a human-following robot using the particle filter.In these applications, data from multiple sensors are combined using fusion techniques to generate a more accurate measurement than from one independent sensor.The Kalman filter is another common fusion algorithm; however, the Kalman filter is only optimal in linear systems with additive, white Gaussian noise.
In general, finding an analytic expression for the estimated states' probability distribution function (PDF) in a nonlinear filtering problem is often impossible.In this case, the EKF may be applied based on the idea of linearizing the system about certain states, then applying a Kalman filter to the linearized model [21].The deficiencies associated with the EKF are: first, significant errors may result from linearization, and it is nontrivial to derive the Jacobian matrices, which may result in implementation difficulties; second, the EKF is not applicable to non-Gaussian noise as described in [21].
The UKF has been a popularized alternative to the EKF, as given by [13], [14], [15].The UKF is also a recursive minimum mean square error estimator just like the EKF; however, the UKF does not require the linearization of a nonlinear motion model.Instead, the UKF uses a set of deterministically selected points propagating through the nonlinear system, known as sigma points, to represent the mean and covariance of the system.Grid-based methods that rely on fixed grid approximation suffer from the curse of dimensionality, as the approximation error depends on the state dimension of the underlying jump Markov linear system [22].Finally, the sequential Monte Carlo method is widely known as particle filtering and relies on approximating the distribution of the states by many weighted samples called particles and is discussed further next.
Particle filters can handle nonlinear, non-Gaussian scenarios, making them ideal for tracking applications.For example, particle filters have been used to fuse data from multiple sensors to track target location over a sensor network that covers an area of interest [23], [24], [25], [26], [27], [28].Furthermore, particle filters can be used in conjunction with inertial sensors for motion tracking applications [10].This application was shown in [25], where fused data from visual sensors and an inertial sensor were used to track the motion of a robotic arm.Sensor fusion is also utilized in navigation applications to provide more accurate positional solutions when GNSS is unavailable.In [26], a hybrid Kalman-particle filter approach was used to produce a navigation solution.While the hybrid approach can provide more accurate results than the purely Kalman-based approach, a complete particle filter solution can provide even more accurate results as the trajectory of a navigating body will no doubt be a nonlinear process.
In general, the particle filtering technique has drawn increasing attention in recent years because of its superior ability to deal with processes with nonlinear models and contain non-Gaussian noise sources [18], [19], [20], [21], [29].Therefore, this article explores the method of fusing outputs from multiple particle filters to obtain more optimal results.The particle filter mathematics are foundational to the proposed algorithm; specifically, the resampling techniques used are derived heavily from the particle filter.Specific details of the particle filter mathematics and the benefits of using the particle filter to interpolate position solutions are discussed in [10].The method proposed in this article is applied to address the high positional accuracy required in SAR applications, as discussed in the next section, with a specific focus of miniaturizing the navigation system.

III. SAR IMAGING
The accuracy of the radar position at each pulse is a critical factor in SAR applications.The quality of the focused image of the target scene is directly related to the ability to accurately record the radar position.Fundamentally, SAR requires the alignment of phase histories returned from the scatterers in the target scene.The alignment has to be coherent on a pulseby-pulse basis, and the phase histories are Doppler processed before coherently integrated to form an image of the scatterers within the illuminated scene.Ideally, GNSS and IMU components are chosen to create a system that have the same, or similar, sampling rates as the pulse repetition frequency (PRF) of the radar.The synchronized sampling rates allow for easy correlation between the radar's transmitted pulses and the position estimates from the GNSS and IMU.In general, a pulsed-Doppler radar architecture is most applicable to this type of system.Pulse synchronization is not as critical for frequency modulated continuous wave (FMCW) radar architectures due to the system's ability to be continuously transmitting and receiving.However, the use of FMCW radar architectures for SAR imaging has its complications such as intra-pulse range migration [30].
Position, navigation, and timing (PNT) solutions utilize a GNSS to provide position and timing and an IMU to provide the navigation.In order to synchronize between the PNT solution and the radar, equally critically important is the PNT data output to be exceptionally accurate.In the event low-quality GNSS and IMU sensors are used, techniques such as autofocus are needed to compensate for the inaccuracies of the PNT system.Autofocus algorithms are designed to calculate the phase errors caused by the SAR system, including the position measurement errors.These algorithms compensate for the errors when forming a SAR image [31], [32], [33], [34], [35].SAR systems that do not use autofocus algorithms perform the phase alignment by relying on the data output from the navigation system [36].As a result, high-accuracy GNSS and IMUs are required for most SAR applications [37].The availability of high-accuracy navigation sensors are typically constrained by the available sampling frequency along with the increased cost, size, weight, and power.The primary consideration is choosing a system that can meet the desired sampling requirements and the secondary constraint is C-SWaP.
Phase plays an important role in SAR image formation as the phase of a SAR waveform is affected by the modulation from the radar transmitter and the two-way free-space wave propagation to and from a scatterer.As a result, the total phase measured by a SAR system is calculated as where φ T (t ) is the modulation caused by the transmitter, λ is the center frequency wavelength, and R(t ) is the range to the scatterer.In most SAR systems, φ T (t ) and λ are designed and known while the range to the scatterer is the parameter being estimated given the estimated relative position of the radar from the PNT system.However, the sensitivity to phase accuracy proportionally increases with the radar carrier frequency and the radar wavelength proportionally decreases.From (1), as the wavelength decreases, the amount of phase change for a given range increases.Thus, any inaccuracies in the estimated range to the target will result in poor phase estimation, which becomes excessively large at high frequencies.Therefore, it is essential to use a high-quality GNSS/IMU to ensure that the relative position of the radar platform is accurately defined such that phase errors are minimized and can be compensated during SAR processing.Ideally, the accuracy required to focus a SAR image requires the errors in range, R(t ), to be limited to λ/16.The phase errors tolerable by the SAR system is estimated as The required range accuracy is limited by the carrier frequency wavelength, and the tolerance for positional errors is increasingly smaller as the carrier wavelength is increases.This is a fundamental constraint for high frequency SAR systems as highly accurate position solutions become increasingly harder to obtain.The desire to increase the frequency of SAR systems is driven by the ability to miniaturize physical apertures and components of the system.However, a SAR system with error-free position solutions can still have issues creating a focused image if the sampling rate of the positions is not synchronized with the SAR PRF.While the above equation allows for the calculation of the total phase change on a single pulse, the radar will transmit and receive several pulses throughout the coherent processing interval (CPI), and thus the range from the radar platform to any given scatterer will change on a pulse-by-pulse basis.Before Doppler processing in traditional SAR algorithms, the individual phase histories, correlated to the range profiles of each scatterer, are aligned along the slow-time radar data such that all the power will compress into a single range bin after Doppler processing.This "alignment" is accomplished by slightly increasing or decreasing the phase change of each fast-time radar data set, which can be precisely determined if the position of the radar platform is known on each pulse of the radar throughout the CPI.This is why it is important for the PNT solution to be accurate but also produce position estimates at the same rate as the radar.If the position estimation data is not available for each radar pulse, the phase for each fast-time radar data set cannot be adequately compensated for and the resulting SAR image will look blurry as the power will smear across multiple bins after Doppler processing instead of being "focused" into a single bin.
As a result, popular interpolation techniques have been introduced in scenarios like this to estimate the position at the radar PRF between GNSS/IMU data samples.However, if the GNSS/IMU data samples are too spread apart temporally, techniques such as the linear or spline interpolation are inadequate in SAR imaging applications.The primary reason these interpolation techniques are insufficient is that they are polynomial interpolations and do not capture the external forces acting on the radar.As a result, the interpolation introduces its own phases errors rather than compensating for them.In order to provide a position solution capable of coherently integrating the phase histories from a SAR scene, an interpolation technique is required that accurately estimates the position of the radar when GNSS/IMU data is not available.Such a technique would allow for pulsed SAR systems to no longer be constrained by the sampling rates of the IMU and GNSS systems.As demonstrated in [10], such a interpolation technique exists and is expanded upon in this work to be applied to multiple IMU systems.

IV. MULTI-SENSOR FUSION TECHNIQUE
The proposed multi-IMU particle filter focuses on fusing the particle clouds from multiple particle filters; each propagated from different sensor measurements.A block diagram of the algorithm is shown in Fig. 1.The proposed method uses sequential importance sampling principles [38], [39], and the performance is evaluated based on the mean squared error (MSE).In addition, the proposed method is compared to the traditional method of combining measurements, such as taking the mean of the particle filter estimations, to compare the performance accuracy of the two different methods.According to the sampling theorem, an increased number of samples allows for a more accurate PDF of the measured variable.As a result, the true value of the measured state can be more accurately estimated.The proposed method uses each sensor measurement as additional samples to more accurately estimate the measured acceleration and, ultimately, the platform's position.The increased accuracy provided by the algorithm allows for smaller and lower-grade IMUs to be used to achieve better accuracy.As a result, the low C-SWaP of lower grade sensors can be leveraged without sacrificing the accuracy provided by higher-grade sensors.
The proposed method is inspired by the SIR principle, where weights are assigned to each particle and then used to resample a set of new particles comprised of more highly weighted particles.The proposed method applies the same methodology by assigning weights to the outcomes of each individual sensor estimation.The sensor weights do not take into account the particle weights computed by the particle filter; thus, the sensor weights are a measure of how accurate the particle filter solution is for any given sensor and are not dependent on the particles used to compute the estimation.However, the sensor weights can be calculated similarly to the particle weights as the principles still apply.
In this article, 1,000 particles are used to resolve the position solution for each sensor.The number of particles used is experimentally chosen to produce the most consistently accurate position estimations for this application.The effects of the number of particles is extensively discussed in literature and more detailed explanations can be found in [19], [20], [40], [41].The weighted estimations are then used as the proposal distribution for a new particle filter, which draws a set of particles from the proposal distribution according to the weighted sensor estimations.A new estimate is then calculated by taking the expected value of the new particle cloud, and the weights are calculated using (3).The new particle cloud collects several samples (N = 50) from the weights, and the respective sensor estimation is used as part of the new particle cloud.By doing so, the sensor with the highest weight is selected more often during the creation of the new particle cloud.
The pseudo-code for the implemented algorithm is illustrated as Algorithm 1.
While the above block diagram and pseudo-code provide an understanding of the fusion algorithm, they do not provide Algorithm 1: Implemented Algorithm.
an intuitive understanding as to how many sensors are needed to achieve a desired positional accuracy.Here, the authors utilize the Cramer Rao Lower Bound (CRLB) to provide a first-order engineering tool to estimate the number of sensors, with known error, to achieve the desired positional accuracy.However, note that the CRLB is commonly used to estimate the minimum variance achievable of the particle cloud in a particle filter as it calculates the lowest amount of variance of an unbiased estimator.In this application, the CRLB does not optimize the unbiased estimator by calculating the lowest amount of positional variance.The idea is to calculate the theoretical error of the fusion algorithm by effectively setting the left side of the CRLB to be the desired positional accuracy.This can be thought of as setting the lower bound of the algorithm to be the desired positional accuracy for the system.The theoretical error contribution from each sensor is calculated from the IMU specifications and used to estimate the number of sensors theoretically needed to achieve the desired position accuracy.Therefore, a derivation of the CRLB is provided.

V. MULTI-SENSOR CRAMER RAO LOWER BOUND APPROXIMATION
The CRLB can be used to approximate the number of sensors needed based on the positional accuracy required.The sensor weights are normalized and used to estimate the effective number of independent samples.The sample space is given by each of the sensor outputs as defined by y 1 , y 2 , y 3 , . .., y s (4) where s is the number of sensor outputs.When using SIR to resample the sensor spaces given a set of weights, the sample space is given by where x is the sample drawn from (4) and N is the number of samples drawn from the sample space S. As a result, the sensor weights are used to calculate the effective number of samples (N s ).The weights correspond to w s (1 ≤ s ≤ N ) and the normalized form becomes The effective independent samples are then calculated using where N s is bound by 1 ≤ N s ≤ N showing that resampling the sensor outputs using sensor weights does not introduce a dependency between sensors.Furthermore, resampling produces independent, identically distributed (i.i.d) samples that are drawn from a discrete density as explained in [20].As a result, the independent samples drawn from the sensor measurement can be applied to the central limit theorem (CLT).
The CLT states that the normalized sum of independent samples can be approximated as a normal probability distribution with a fixed mean and variance.This assumption allows for the use of the simplified Gaussian form of the CRLB [42].
As previously mentioned, the CRLB establishes the lowest amount of variance of an unbiased estimator; however, for the proposed algorithm, the CRLB will express the lowest amount of variance that is expected in the navigation solution when fusing multiple sensors [43].In general, the CRLB is calculated using where the denominator is the Fisher information matrix.
The sensors in a navigation platform produce independent samples of the platform position and the SIR method uses the principle of weighted sampling described by (7), where the solution from each sensor is approximated as where x is the true position of the platform and w[n] is white, Gaussian noise with a given variance σ 2 .Given the sensors are independent, identically Gaussian distributed with the expected value being the true platform position, the probability distribution for the sensor estimates can be expressed by where N is the number of sensor solutions.The first step in evaluating the Fisher information matrix is taking the first derivative of the natural log of (10) shown as Evaluating the summation and taking the second derivative simplifies the solution as follows Combining the solution obtained in (V) with the definition of the CRLB expressed in (8) yields (13) establishes the variances of the positional estimates from each sensor, which is bound by a single sensors positional variance divided by the number of sensors used.This is applicable to a single time step and can be used for a back of the envelope calculation for the number of sensors needed.Note, this simplification assumes the variance of each sensor is identical.It is important to note that this technique provides an estimation of the number of sensors an application may require.The method of using this CRLB estimation is based on the approximation that the sensor errors can be modeled as white Gaussian, which is a general approximation and does not fully model the true effects of the sensors chosen.A simulation of 1,000 measurements for a single time step is used to evaluate the CRLB technique.These measurements are resampled using the SIR technique with a variable number of sensors from 1 to an maximum of 100.The sample variance of the position solution is calculated across the sensors.The results are compared against the sample variance of a single sensor divided by the number of sensors used in the SIR technique and is shown in Fig. 2. As a result, the plot shows excellent agreement between the results of the SIR technique and the theoretical estimation using (13).
Additionally, the MSE of the SIR technique as compared to GNSS is shown in Fig. 3.It can be seen that the addition of more sensors results in a decrease in MSE compared to the error of the GNSS measurement used to calculate the sensor weights.As additional sensors are added, the solution converges to the GNSS measurement.As the particle filter predicts and resamples, the MSE will continue to decrease.In Fig. 3, the results are for a single time step, and it is illustrated that the accuracy for a single timestep can be significantly increased using additional sensors and the proposed SIR technique.
In Fig. 3, the error from the fusion technique is improved as sensors are added; however, the MSE does not converge to GNSS completely because the MSE is only shown for a single time step.As a result, the sequential benefits of the SIR technique are not displayed.Moreover, it should be noted that the law of diminishing returns is likely to contribute when designing multi-IMU navigational systems using the proposed method, as illustrated in Figs. 2 and 3. Thus, there is a point where the additional sensors do not provide a significant accuracy increase to the final position estimated.In general, adding additional sensors at some point will become unfeasible due to other constraints in the desired application.Furthermore, the simulated results are selected to simulate out to 100 sensors to fully observe the convergence of the CRLB approximation.In practice, this number is very likely unfeasible and the number of sensors utilized in instrumented systems will likely be much less than 100.

VI. NON-GAUSSIAN NOISE SENSOR MODEL
In order to calculate the CRLB of the sequential Monte Carlo, a distribution for the IMU sensors needs to be modified from (9) to more accurately represent sensors with non-Gaussian distributions.Assuming an expected value and variance can model the distribution, the variance of the positional solution needs to be estimated.The expected value can be assumed to be the true position of the platform.The variance of the positional solution can be viewed as the errors in the position estimates from the IMUs.In this case, the positional errors are caused by integrating the errors from the accelerometers (a e (τ )), the initialization error in velocity (v e (0)), and the initial position error (P e (0)), calculated by v e (0)dτ + P e (0).(14) The errors in the acceleration measurement from the IMUs can be calculated by accounting for the errors in measured acceleration, the errors in rotating the measured accelerations to the right coordinate frame, and the errors in the gravity and earth rotation models using In (15) where θ x,e , θ y,e , θ z,e are the estimated x, y, and z errors of the platform's attitude, respectively.The IMU's error effects can be estimated from the specifications provided by the sensor's manufacturer.In addition, bias instability (b), white noise (N), and random walk (RW ( f ) = K 2 /(2π f ) 2 ) are also contributing factors to the IMU errors.Moreover, the scale factor errors at a given acceleration (B(a t )) are combined with the misalignment (θ m ) and non-orthogonality (θ n ).All these variables are used to calculate the acceleration errors using (17) Lastly, the IMU errors can be used with the true acceleration (a t ) and the duration (T ) during which the acceleration is experienced to calculate the positional standard deviation calculated using

VII. CRLB WITHOUT GAUSSIAN ASSUMPTION
For the SIR method, the CRLB is insufficient in estimating the variance because the method is recursively used to represent the posterior distribution.As a result, the CRLB for a single time step only provides a first-order approximation of the positional accuracy of the SIR method.The variance of the SIR method can be more accurately approximated by using the sequential Monte Carlo approach.The variance of the Monte Carlo solution can be estimated using (19) In (19), the probability density of the sensor estimates is represented by π n ( xk ), where n is the sensor and k is the time step.Furthermore, the proposal distribution is represented as q k ( Xk ) and is the distribution used to calculate the weights of each sample during the SIR process.Equation ( 19) captures the recursive nature of the SIR method where the first time-step, This distinction allows for the distribution of the previous time step to be utilized as the proposal distribution in the next time step.The proposal distribution for time steps two through N are described by A full derivation of the sequential Monte Carlo variance is presented in [40].In the multi-IMU fusion technique, q k is used to calculate the sensor weights and is a distribution representing the GNSS sensor.Furthermore, π n ( xk ) is a distribution representing the positional solutions from the IMUs.As a result, solving (19) for N gives an estimate of the number of sensors needed to accomplish a desired positional accuracy Var{ X } and the desired number of time steps k.Using ( 13) and ( 18), the variance of the fusion solution can be calculated to have a lower bound of σ 2 /N.As a result, the effects of the additional sensors can be estimated.Given a variance in the positional accuracy of all the sensors used, the number N can be varied to achieve the desired positional accuracy needed for the desired application.This result can be used as a design tool to understand and select the appropriate sensors to satisfy the C-SWaP and accuracy constraints in estimating the position of a navigational platform.
Using the CRLB to estimate the number of sensors required to achieve a specific accuracy also allows for an estimation of the multi-IMU navigation systems C-SWaP.As a result, the estimations allow designers to quickly determine if the chosen sensors are beneficial in decreasing C-SWaP while obtaining the desired positional accuracy.Other considerations also apply to the design of multi-IMU navigation systems, such as the computational bandwidth and data storage associated with the additional data from the added sensors.This obviously will play a significant factor in the number of sensors that can be effectively added.In addition, the computational time to achieve a positional solution may be unfeasible for certain applications.Finally, space constraints may also be another factor in the design of a multi-IMU system, as the constraint for lower C-SWaP may be based on limitations in the platform size; especially when considering small airborne platforms.

VIII. SIMULATION
In order to emulate a realistic navigation scenario, a flight path, shown in Fig. 4, is simulated using three IMUs integrated with a GNSS solution.The IMU parameters used to simulate this scenario are derived from the Analog Devices ADIS16465 precision IMU and are specified in Table 1.For this simulation, each particle filter was initialized to the starting position and velocity of the navigating body, with the variance of the particle cloud equal to the uncertainties of the Analog Devices IMUs.These uncertainties are derived from the IMU specifications as detailed by [37]Next, each IMU measurement is transformed into the local navigation frame (North, East, and Down) using measured gyro angular rotation rates.In addition, acceleration values from the IMUs are used to resolve the body's velocity in the local navigation frame.The velocity is then used to update the position in geodetic coordinates (latitude, longitude, and altitude).Then the positions are resolved using a particle filter fused with simulated GNSS measurements.Finally, the error between the algorithm's predicted position and the actual simulated path is calculated for each sensor independently and is shown in Fig. 5 as a function of time.Note that the three IMUs used in this simulation had the same noise and biasing characteristics but yielded different error values.This result indicates that the inaccuracies of an IMU are not deterministic, even if the IMU specs are the same.Therefore, the individual IMUs can be treated as independent samples, which the fusion algorithm requires.
In simulation, two different flight paths are defined, and the individual sensor position estimation capabilities are investigated.Here, the individual sensor results are fused using the proposed algorithm for each simulated path.The error is calculated between the simulated "truth" position and the estimated position from the proposed fusion algorithm.Furthermore, the same error calculation is done with the position estimations from a traditional fusion method of taking a normal average of the multiple sensor estimations.These calculations allow not only a comparison of the enhanced positional estimation between the proposed algorithm and individual sensor estimations, but also a comparison between the proposed algorithm and a traditional averaging fusion method.
The error of the traditional and proposed fusion methods are shown in Fig. 6 for the first simulated path and associated sensors.Note that traditional fusion substantially reduces positional error compared to Fig. 5; however, the proposed fusion method is approximately three times more accurate than the traditional method.It is important to note the traditional method here is simply taking the average of the multiple sensors.Furthermore, the proposed method of resampling the outputs of the sensors instead of taking the average significantly improves the accuracy of the positions estimated.It is this finding that allows the method proposed in this article to be capable of producing the position accuracy required for SAR.The same comparison is shown in Fig. 7 for the second simulated flight path and the associated Analog Devices IMUs.Again, there is a significant reduction in error using a traditional fusion method (almost five times reduction), but the proposed fusion algorithm provides approximately an additional two times improvement.Based on the results in Figs. 6 and 7, the proposed fusion algorithm provides excellent position estimations of a navigating body compared to individual sensor capabilities and traditional averaging fusion algorithms.
Given the proposed algorithm's excellent ability to estimate position in simulation, the final task is to compare the performance between a single navigational grade IMU and the fused output of the three precision MEMS grade Analog Devices IMU modules.The navigational grade IMU chosen for this  2. Due to the increased C-SWaP of navigational grade IMUs, their usability in SWaP-constrained environments, such as airborne platforms, is drastically limited.Therefore, the position, navigation, and timing (PNT) module could be greatly limited in these situations to smaller, lower-grade IMU modules.However, the position estimation accuracy of the system will more than likely need to be maintained, albeit these constraints and the use of multiple reduced C-SWaP and quality IMUs with the proposed fusion algorithm could be an ideal solution.
The second flight path is simulated with the single navigational grade NovAtel IMU and a GNSS solution, and the data from each are fused with a particle filter to yield the estimated position.Next, the positional error between the simulated "truth" path and the estimated path by the fused GNSS/navigational grade IMU is calculated as a function of time and shown in Fig. 8.The error for a single ADIS16465 and three ADIS16465 IMUs are shown as well.It is apparent that the NovAtel IMU has a lower error when compared to a single precision MEMS grade Analog Devices IMU as expected.However, the three ADIS16465 IMU fused solution resulted in an error on par with the NovAtel unit.Based on the IMU specifications shown in Tables 2 and 1 as inputs into the methods mentioned in section V, three sensors is determined to be able to provide similar positional quality compared to the navigational grade NovAtel IMU.This result indicates that a highly-precise PNT solution can be achieved by fusing multiple lower-grade IMUs using the proposed fusion algorithm.

IX. INSTRUMENTATION
In order to validate the simulated results, an experimental setup is constructed consisting of the NovAtel system and three ADIS 16465 IMUs attached to a vehicle that is driven through an urban environment as seen in Fig. 9.
The instrumented setup is shown in Fig. 10 where the No-vAtel system is the two grey boxes comprised of the IMU and the ProPak6 TM receiver.In contrast, the Multi-IMU system is the gold box and the three small metal boxes on top of the NovAtel IMU.The gold box is a Remote Sensing Solutions (RSS) control and timing unit (CTU) which acts as a receiver for the three ADIS16465 IMUs.The three IMUs are attached to their respective evaluation boards and communicate with the CTU by ribbon cables.The CTU and the three IMUs are small enough to be deployed on top of the Novatel IMU and weighs significantly less than the NovAtel IMU alone.In comparison, the Multi-IMU system, with the CTU receiver, is a reduction in C-SWaP with respect to the NovAtel IMU alone.The C-SWaP reduction is further increased when compared to the NovAtel IMU and the NovAtel receiver as a system.While the biasing and error characteristics of the No-vAtel unit are significantly smaller, the unit has an increased order of magnitude in cost, size, and power consumption and is three and a half times heavier compared to the three Analog Devices solution.Therefore, it is clear that there is a C-SWaP vs. performance trade-off between the different PNT solutions.
The instrumented system collects GNSS latitude, longitude, and altitude from the NovAtel system.In addition, the angular velocities and acceleration values from the NovAtel IMU-ISA-100 C are collected from the NovAtel synchronized position attitude navigation (SPAN) receiver.The angular velocities and acceleration values from each of the three IMUs are collected using the RSS CTU.The data collected from each sensor are timestamped by the CTU and the NovAtel SPAN receiver.In addition, the experiment also utilizes an RTK GNSS with centimeter-level accuracy to capture the path the vehicle traverses.The RTK measurements are used as a highly accurate third-party ground truth of the path traversed to measure each IMU configuration's accuracy.
The data from the NovAtel IMU is fused with the NovAtel GNSS data using the traditional particle filter fusion method.For a fair comparison, the three ADIS IMUs are also fused with the NovAtel GNSS measurements using the proposed fusion method.Fusing each IMU configuration with the same GNSS measurements allows for the uncertainty of the GNSS measurements to affect all the IMU configurations equally.As a result, the differences in errors from each IMU configuration are primarily sourced by the errors from the IMU and the fusion method.As seen in Fig. 8, the simulations show that three ADIS 16465 IMUs can provide comparable performance to the NovAtel IMU.Therefore, the experiment is designed to confirm these results from the theory and simulation.

A. CALIBRATION AND ALIGNMENT
The instrumented setup is calibrated using the alignment method from the NovAtel system software.The instrumented results are powered on in advanced to allow for the alignment process to complete.After the system is aligned, the system remains powered for the duration of multiple SAR data collections.The idea is to preserve the calibration for the navigation duration rather than performing a calibration before each SAR data collection.In post-processing, the time stamps from the NovAtel system, the multi-IMUs, and the SAR pulse timing is interpolated and synchronized using the method in [10].The initial measurements from the multi-IMU system are then compared with the initial measurements from the NovAtel to achieve calibration synchronization between the two systems.As the multi-IMUs are mechanically affixed to the NovAtel system, the synchronization of the two initial measurements provides the same alignment for both systems.As the target application is sensitive to accumulated errors during the SAR collection, the synchronization of the calibrations is sufficient.

B. MEASURED RESULTS
The results of the experiment for the NovAtel IMU, single ADIS 16465 IMU, and the three ADIS 16465 IMUs are shown in Fig. 11.The results are compared to the RTK GNSS, and the error of each IMU configuration is calculated.Fig. 11 shows that the three IMU configuration performs the best, followed by the NovAtel and the single ADIS configuration, respectively.These results are expected as the NovAtel system is a more accurate sensor compared to a single ADIS 16465.Similarly, the NovAtel system and the ADIS start with similar accuracy, but the single ADIS solution begins to "walk" 20 seconds faster than the higher quality NovAtel.Therefore, the three IMU configuration with the proposed fusion method improves the initial accuracy of the navigational solution and minimizes the walking effects of the sensors.

X. SAR MOTION COMPENSATION DEMONSTRATION
To demonstrate the proposed technique's applicability to motion compensation in SAR imaging, the instrumented  multi-IMU navigation system is used in conjunction with the NovAtel navigation system to perform motion compensation on data captured using a Ku-band (15 GHz) SAR system.The SAR data is captured from a radar system placed on a moving cart with both the instrumented multi-IMU system and the NovAtel system.
The radar system used in the test is a custom-built system designed at the University of Oklahoma and is a modified version of the radar system described in [44].The system operates at a carrier frequency of 15 GHz, giving a carrier wavelength of 2 cm.Based on the heuristic criteria in (2), the required motion compensation accuracy to produce a focused SAR image is 1.25 mm.The radar's pulsed waveform is a linear frequency modulated (LFM) pulse with 100 MHz of bandwidth, giving a radial range resolution of 1.5 m after the pulse compression operation.To capture the radar data and provide pulse synchronization to the NovAtel unit, a Remote Sensing Solutions (RSS) ARENA413 software-defined radio (SDR) is used.
The scene for the experiment comprises two vehicles acting as reflectors, as well as some foliage.The trucks are placed at 130 m and 145 m in range from the radar.An image of the SAR scene is given in Fig. 12.The SAR system is moved on a cart over a synthetic aperture length of 2 m such that the azimuth resolution at 140 m is approximately equivalent to the range resolution, giving roughly square resolution cells in the resulting SAR images.Both the NovAtel and the instrumented multi-IMU systems are placed on the cart alongside the radar to capture its motion profile during the radar capture.Both systems utilize the same GNSS subsystem to again ensure any GNSS measurement uncertainty is applied to both systems fairly.
A separate SAR image is formed for each system using the same radar data but with motion compensation performed based on the outputs of each separate navigation system.The images are formed using the backprojection algorithm.The SAR images formed with each navigation system are shown in Fig. 13.
Visual inspection shows that the multi-IMU system can better focus the image, leading to a lower sidelobe level and a more ideal sinc-like point spread function in the azimuth dimension.The tighter point spread function is a result of less severe time-varying motion compensation errors.The results obtained by both systems show that the multi-IMU is capable of producing high accuracy position solutions on par with the high-end NovAtel IMU.In addition, the system is capable of producing PNT for SAR applications when paired with the interpolation method.The results further show that the fusion of multiple lower grade IMUs can be used as a SAR PNT solution.
The SAR images acquired using the multi-IMU technique successfully address the challenge of INS inability to provide accurate position updates at the PRF of SAR.Additionally, the multi-IMU technique successfully improves the position quality of the INS while simultaneously minimizes C-SWAP.The capable improvement of the multi-IMU method is demonstrated by the focused SAR images acquired from the ground test where the multi-IMU technique is outperforming the high-grade NovAtel INS.As a result, the findings of this work support the ability to miniaturize SAR systems through the ability to decrease C-SWAP of the required INS integration in SAR systems.

XI. CONCLUSION
A significant limitation in miniaturizing navigation platforms is the issue of requiring a single large C-SWaP IMU to produce an accurate result.This article proposes a particle filter approach to fusing multiple, co-located, low C-SWaP MEMs grade IMUs designed explicitly for lower C-SWaP, nonlinear, non-Gaussian scenarios.The results show that the proposed method is capable of providing similar or better accuracy than the traditional single high accuracy IMU.As such, the presented CRLB enables the instrumentation of the proposed method by allowing engineers to design and estimate the number of sensors needed to achieve the desired positional accuracy.The CRLB presented allows for engineers to connect the data sheet specifications to the desired positional accuracy needed for the designed application.The simulated results confirm that the SIR technique with multiple lower-grade IMUs produces similar accuracy measurements as a single large C-SWaP navigation grade IMU module.
The SIR method is dependent on the fact that multiple MEMs grade IMUs produce statistically independent samples that can be resampled much like the particles in a particle filter.Furthermore, the SIR method is proven through instrumentation and measurement where three MEMS grade IMUs are compared to the navigational grade NovAtel IMU.In addition, an independent RTK GNSS system is used as a highly accurate "ground truth" to estimate the absolute error of each navigational solution.The instrumented results illustrate that three ADIS16465 IMUs can achieve similar positional accuracy of the NovAtel IMU, thus, confirming the theory, algorithm, and simulations presented.The authors acknowledge the particle filter is computationally complex and can be a limiting factor for some real-time applications; however, the proposed method was designed with post-processing in mind.
Finally, the images formed from the instrumented SAR systems directly show that the multi-IMU navigation solutions provides a more focused image than the single-IMU NovAtel system.The improvement of using the multi-IMU system is observed by the point spread function showing a more sinclike point spread function in the azimuth dimension of the focused SAR image.The more ideal point spread function results in a more focused image of the target scene.As a result, the instrumented multi-IMU SAR system further validates the methods and theory proposed.Lastly, the instrumented system used in this article results in a combined three orders of magnitude reduction in C-SWaP and enables the miniaturization and instrumentation of future low C-SWaP navigation solutions.

FIGURE 1 .
FIGURE 1. Block diagram of the proposed algorithm.

FIGURE 2 .FIGURE 3 .
FIGURE 2. The variance of the multi-IMU fusion method compared to the simulated sample variance for a single time step.

FIGURE 4 .
FIGURE 4. Simulated path for the three IMU scenario.The path is generated using Acienna's gnss-ins-sim python library.This path simulates a more complex path that has changes in latitude, longitude, and altitude.The python library simulates the gyroscope measurements as well as the accelerometer measurements to fully simulate measurements from an IMU.

FIGURE 5 .
FIGURE 5.The calculated error of three ADIS16465 sensor simulations compared to the true path.

FIGURE 6 .
FIGURE 6. Error of both algorithms.The blue line is the traditional method of taking the mean and the red line shows the results for the proposed method.

FIGURE 7 .
FIGURE 7. Error of three IMU scenario.The blue line represents the traditional method of taking the mean and the red line shows the results for the proposed method.TABLE 2. Novatel IMU-ISA-100 C Specifications

FIGURE 8 .
FIGURE 8. Error of NovAtel unit compared to the ADIS16465 IMUs.

FIGURE 9 .
FIGURE 9.The path traversed by the deployed NovAtel and ADIS system.This image is a screenshot from Google Earth.

FIGURE 10 .
FIGURE 10.Instrumented setup of the NovAtel system and the Multi-IMU system.

FIGURE 11 .
FIGURE 11.Error of the positional accuracy in meters for the different measurement configurations.

FIGURE 12 .
FIGURE 12. Image of the radar scene.The truck targets are highlighted in red and the trees in the scene are highlighted in green.

FIGURE 13 .
FIGURE 13.A SAR image formed using the NovAtel system (left) and the multi-IMU system (right) for motion compensation.Both images are plotted in dB scale.