Tightly Coupled Integration Design for a Magneto-Inertial Navigation System Using Concurrent AC Magnetic Measurements

Recently, magnetic measurement-based pose estimation technology increasingly attracts the researcher’s attention as a complementary navigation technology. The main purpose of introducing magnetic fields is to secure a reliable pose estimation in non-line-of-sight (NLOS) environments where satellite navigation is unavailable or vision, lidar, and other auxiliary sensors that make use of signals of opportunity cannot provide robust performance. Unlike the previous researches, this paper proposes a new method for implementing the integrated magneto-inertial navigation system based on three-dimensional AC magnetic field vectors from concurrent measurements. Instead of taking an independent pose estimation result from the magnetic positioning system, the proposed algorithm investigates a direct integration scheme that employs the raw magnetic field vectors from concurrent measurement sources. A theoretic formulation is developed, which enables to establish the integrating framework with the INS mechanization. The presented algorithm is validated through trajectory experiments in a small-scale test platform. As a result, it is demonstrated that the computational efficiency is substantially enhanced to enable real-time implementation of the proposed magneto-inertial navigation system, while the estimation performance has sub-centimeter accuracy in a 3D dynamic experiment.


I. INTRODUCTION
The positioning technology is one of the main issues in the navigation field, and there have been many related studies. As a representative technology, GNSS is utilized in many industrial and scientific fields. GNSS system that is based on electromagnetic radio wave signal, however, cannot guarantee its accuracy performance under weak GNSS environments such as indoor, underwater, or near a building. Many alternative methods have been developed to resolve this limitation. One of the related studies includes optical measurements such as vision sensor, laser range finder, or lidar that generates aided information for positioning. [1]- [3] However, despite the abundance of available data in such optically benign environments and the variety of developed image processing techniques, limited application is allowed to environments where the acquired image data is very susceptible to light condition or the line-of-sight (LOS) The associate editor coordinating the review of this manuscript and approving it for publication was Di He .
is not guaranteed. WLAN or UWB is also widely used for indoor navigations with a non-LOS condition. According to their typical electromagnetic characteristics, the relative distances between the signal emitters and receivers are measured by RSS (received signal strength) or ToF (time of flight) [4], [5]. Then a multilateration method is usually applied to obtain position fix via a typical least square or nonlinear estimator. A different approach based on RSS is to match the received strength data with the pre-obtained signal footprint. As is typical in the matching-based positioning scheme, predicting relative distance is not needed, and determining position is free from requiring a strict LOS condition. On the other hand, positioning performance highly depends on a prerequisite map accuracy, whereas the computational efficiency of the matching algorithm itself may degrade due to map complexity and data processing burden.
As an alternative positioning strategy in an NLOS environment, the localization method using magnetic measurement has attracted the researcher's interest in recent related studies [6]- [16]. At first, the perturbation map of the VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ static geomagnetic field is utilized for positioning research [7]- [9], and a moving magnet position is detected by sensor array configuration [10], [11]. These methods fundamentally measure DC magnetic field that has advantages in practical system implementation. Typically, DC magnetic field can be generated in a narrow space by using a strong permanent magnet or pulsed electromagnet. However, owing to its field coverage characteristics, DC magnetic field-based positioning scheme has been limitedly applied to the medical or rehabilitation cases. Moreover, due to the rapid field intensity attenuation, magnetic field distortion potentially appears near static magnetic impediments.
To improve the aforementioned shortcomings of the DC-based technique, an AC magnetic field-based positioning concept is investigated alternatively. Compared with the DC magnetic implementation, the AC magnetic field can propagate further even with the same transmitting power, thus provide an extended measurement coverage. Also, as each AC magnetic field can be uniquely extracted from the concurrent magnetic field signals with multiple carrier frequencies, an enhanced estimation performance can be achieved with the help of the augmented measurement dimension.
Positioning methods employing time-varying magnetic fields have been reported increasingly in recent studies [12]- [16]. Among these, some study cases employ TDMA (time division multiple access) or FDMA (Frequency division multiple access) approaches for signal extraction from concurrent magnetic fields [12]- [14]. In contrast, the AC modulated signal is processed to provide 3D pose determination using the magnetic measurement from a three-axis magnetometer or multiple receiver coils [15], [16]. Furthermore, as a distinguished sensor fusion for navigation purposes, magnetic field measurement is further integrated with other navigation sensors such as INS (inertial navigation system) [17]- [19], [21] and GNSS (Global Navigation Stellate System) [20]. As mentioned earlier, integrated navigation using DC magnetic measurements revealed to yield low accuracy in correcting pose estimation. Besides, the drawback of limited measurement range or high sensitivity to the ambient magnetic noises around the test environment was demonstrated. In contrast, the author's previous study introduced an innovative integrated navigation scheme combining inertial sensors and magnetic positioning system [21]. In this study, the proposed MPS(Magnetic Positioning System) first computes the position and attitude estimation through developing a grid map representation model, and then the predicted pose was combined with INS through the integration filter framework. Despite positioning error performance of sub-centimeter level, loosely coupled integration scheme and burdensome magnetic pose determination algorithm has been proved as a significant obstacle to the real-time implementation of the system.
In this background, this paper proposes a new integration mechanism that directly combines inertial and AC magnetic measurements for the pose estimation. For this, a tightly coupled sensor integration scheme is developed for a sequential filtering process, which enables the algorithm's real-time implementation. The performance of the proposed algorithm is analyzed through a dynamic motion test with 2D and 3D trajectory. From experimental results, the advantage of the tightly coupled integration scheme is verified by demonstrating the real-time performance of the integrated navigation solution.
The main contributions of this study are summarized below. First, a new spatial representation model based on concurrent AC magnetic fields was developed, which enables to formulate pose-dependent magnetic measurement. This suggests a fundamental difference from the conventional multilateration positioning technique based on the magnetic field attenuation effect. Also, this paper exploited a novel integration scheme that incorporates the inertial measurement to secure the accuracy and continuity of the magnetic field-based navigation. In particular, by deriving a filter formulation that combines AC magnetic measurements directly, the computational inefficiency of the loosely-coupled integration approach in previous research is overcome. Finally, the presented test results of the magneto-inertial navigation contributed to demonstrate the real-time performance of the proposed method.

II. AC MAGNETIC FIELD BASED MAP DATA
In this section, system configuration and spatial representation model for constructing AC magnetic map is presented. Figure 1 shows the conceptual diagram of the proposed system configuration. As observed in the figure, electromagnets excited by different carrier frequencies are spatially deployed for constructing fundamental AC magnetic grid map. Each AC magnetic field is sensed by a three-axis magnetometer, which provides the magnetic field signal in analog voltage with the superposed magnetic signal waveforms. For a fast data processing, an embedded DSP processor is used for analog-to-digital converting with fast sampling rate higher than 1kHz. In our practical implementation, four electromagnets with carrier frequencies of f 1 (40 Hz), f 2 (50 Hz), f 3 (60 Hz), and f 4 (70 Hz) are used for instance. Since the magnetometer has three-axis, these raw signals have three components, (x magn , y magn , z magn ). Equation (1) represents a model of the raw signal sensed by the x-axis of the magnetometer, representatively. The subscript 'b' at (1) means sensor output is obtained from the B-frame (i.e., body frame).

A. MAGNETIC VECTOR EXTRACTION
where, f i means each carrier frequency, m f i means magnitude of magnetic field intensity for each carrier, ω i = 2 · π · f i , and m 0 represents a DC component of magnetic signal. To extract each AC field vectors from the raw signal, the reference sinusoidal signal is employed, which has the same frequency as the exciting AC signal of each electromagnet for generating magnetic field.
Equation (2) depicts a component of reference sinusoidal signal with frequency of ω 1 .θ 1 represents the locked phase for demodulation.
Equation (3) shows the integral result of the multiplication between magnetic signal and reference sinusoidal signal. As is similar with a typical demodulation procedure, by taking an integral window with the matched time interval and low pass filtering, the informative magnitude of magnetic field, m f 1 ,b can be acquired. Therefore, the magnetic field, m f 1 ,b , can be obtained as the following equation by removing all the harmonic terms of f 1 .
where t represents integration interval. In (4), integration interval is assigned as the inverse of the maximum common factor of all carrier frequencies. And, θ 1 andθ 1 are the AC magnetic signal phase and its detected phase for demodulation, respectively. If the phase condition is matched (i.e., θ 1 =θ 1 ), the result of the integral is maximum. The integral is computed numerically through digital processing, where phase sweeping method for f 1 is employed for best approximating the phase estimate. Note that phase resolution during sweeping process may cause a magnetic field magnitude uncertainty, while excessively small phase resolution results in computing burden. Note that all the frequency terms except ω 1 and bias, m 0,b in (1) are removed through demodulation process. This feature serves as an essential advantage of adopting an AC magnetic field, which enables to neglect hard iron effect compensation. Soft iron effect due to sensor's inherent scale factor nonlinearity may still remain as error source. A preliminary calibration process to compensate for soft iron effect is adopted in this study, yet its detailed description is out of scope and omitted here. Therefore, without loss of generality, a 3D magnetic field vector can be repeatedly obtained using the demodulation and calibration process for each carrier frequency. Equation (5) shows the magnetic vector representation composed of 3 dimensional elements, which is denoted by M i .
where m i = 2 t · t+ t t x magn · r i dt, n i = 2 t · t+ t t y magn · r i dt, and l i = 2 t · t+ t t z magn · r i dt. Note that subscript i denotes each carrier frequency.  Figure 2 illustrates overall conceptual block diagram of signal demodulation from multiple AC magnetic measurements. As shown in Fig. 1, four AC magnetic carrier signals are concurrently superposed in the same space. For obtaining each magnetic field with different carriers, raw measurement VOLUME 8, 2020 is demodulated through numerical phase search using its corresponding reference signal, r i . Consequently, four magnetic vectors can be obtained as M i , i = 1,2,3,4, from the raw AC magnetic measurements. After constructing the grid map with concurrent magnetic vectors, a 3-dimensional magnetic field model at the position of interest is investigated through the following interpolation strategy.

B. MAP REPRESENTATION MODEL VIA INTERPOLATION
In this subsection, map representation model via interpolation is introduced. It is notable that a similar map representation model is employed in author's previous work [21], which computes standalone pose estimation via optimal map matching approach. This paper, on contrast, suggests incorporating AC magnetic measurements into the integration filter framework by exploiting the spatial representation model of grid map.
In representing the AC magnetic field, the interpolation model is fundamentally adapted to cover magnetic data at an arbitrary position in the grid map. The proposed scheme enhances the effectiveness of map acquisition process with a proper grid distribution, while interpolation model successfully provides geometric formulation of magnetic field vector in the constructed map.   Now the spatial representation model is formulated as below. First, consider the 27 grid points in the map surrounding the objective position (uvw). Figure 4 shows the relation between the objective position and the grid map, where the center grid is located at position (222). Initial center grid is determined such that it minimizes the norm of vector difference between magnetic measurement and grid map data. Also, while each grid point is composed of four 3D magnetic vectors with each carrier frequency, the same computing process can be repeated using the general formulation as shown below.
Without loss of generality, let m(uvw) denote the magnetic measurement in the X-axis at a position (u, v, w). Computing magnetic interpolation in terms of position is suggested in the order of coordinate axis as follows. First, in the Z-axis, it is described as Next, in the Y-axis, it is described as follows; Then in the X -axis, the computed magnetic vector is related with Note that D is a constant matrix since the distance between the map points is initially determined during grid map generation. Thus, based on a second order interpolation model, the computing process of (6)-(8) yields the predicted magnetic field model at position (uvw) in the X-axis. Similarly, the same operation on the Y-axis and Z-axis generates n (uvw) and l (uvw), respectively. Finally, a repeated calculation with the interpolation model can be applied for each carrier frequency. The resulting 3-dimensional magnetic field representation model for each carrier frequency is generally described as with an arbitrary position, (uvw) in the grid map. Note that analytic magnetic vector from grid map is denoted by, which is distinguishable from the true measured vector, m i . Also since the map data is constructed by the magnetometer aligned with a local horizontal plane, all the magnetic vectors in grid are aligned with the N-frame (i.e., navigation frame). Consequently, the spatial representation model in (9) is represented in the N-frame.
In author's previous study, the measured AC magnetic field is iteratively compared with the magnetic grid map data to satisfy a predefined threshold with the help of magnetic interpolation model. Then the estimated pose information via magnetic data analysis is further integrated with the INS in a loosely coupled strategy. Despite a sub-centimeter error performance in positioning, however, excessive computing resource is required for searching candidate grids and resizing grid to terminate iterative matching process. This resulted in difficulty for achieving realtime performance of the previous pose estimation method.
To overcome the raised performance limitation, this paper suggests a new estimation filter scheme employing a tightly coupled integration with the magnetic measurements. Compared with the iterative matching method in the previous work, the proposed method takes advantage of analytic linearization on the spatial representation model of magnetic vectors. Thus, measurement update in filter is directly implemented using the derived formulation instead of operating burdensome numerical search in the candidate space. Detailed description is presented in the next section.

III. INS AND MAGNETIC MEASUREMENT INTEGRATION
In this section, mathematical development for system integration is presented. For this, INS model with modified state variables is introduced, which is further developed to construct efficiently filter formulation. Especially, an approach with the mixed attitude representation is adapted by taking both quaternion and 3-component vector form to avoid singularity and to secure computational efficiency.

A. INS MODEL AND LINEARIZATION
In this subsection, INS error model and state propagation through IMU output are developed. The adopted INS model includes navigation states of position, velocity, attitude and bias of accelerometer and gyroscope as below, where P, V , represent position, velocity and attitude error in the N-frame, respectively. b acc , b gyro represent sensor biases in the B-frame. Note that navigation state is composed of 15 elements and further developed for linearized INS model. Note that attitude error in state model takes three angles form with small angle assumption for a simpler development of filter formulation. However, in practical filter implementation, a quaternion model is adapted to update attitude with respect to reference navigation frame.
Equation (11) introduces the attitude error in the threeangle form. φ and (α, β, γ ) are a scalar value and the direction of the attitude vector, respectively. Intuitively, three component vector implies the attitude error from the estimated attitude to the true attitude. In real implementation, the attitude in quaternion form is parallelly processed using the computed three angle components. Thus, attitude quaternion can be propagated and updated at every epoch, and used to provide the DCM (direction cosine matrix), while the attitude error model depicted by the special orthogonal group in three dimensional space(i.e., SO (3)) is used effectively for implementing the linearized filter structure. The parallel usage of dual attitude representation has advantages of nonsingular parametrization while retaining a reduced dimensionality in filtering application [22,24].
First, consider the propagation model of velocity and position in INS model as below.
where, a n = C n b · A b − b acc,k − g n and g n = 0 0 −g T .The output of the accelerometer A b is calibrated by the accelerometer bias b acc,k , and the DCM from b-frame to n-frame C n b is generated by the function f DCM q k . g n is a normal force by the earth gravity. Considering the IMU error property, constant bias is assumed for accelerometer and VOLUME 8, 2020 Now attitude model relating the attitude error in the form of SO(3) with a quaternion computation is described. First, consider an attitude quaternion form, with four elements q 0 q 1 q 2 q 3 T . The first element of the quaternion q 0 is a scalar term and the other three elements q 1∼3 representing vector terms.
where q − k+1 is a priori update of quaternion at k + 1, q ω = f qua p q r T · t , and p q r (14) shows the quaternion attitude propagation, whereas the output of the gyroscope ω b is calibrated by the gyroscope bias b gyro,k . t equals to IMU update period. f qua (·) converts the three-component angles to quaternion, and means the quaternion multiplication. The geometric relation of the 3D angle is shown at (15).
In (15), * means a quaternion conjugate,ˆmeans the state estimate, and ξ is the attitude error expressed in quaternion from. Then a posteriori update is computed by (16) where, qb In a typical INS model, the dominant error source of attitude mostly comes from gyroscope bias. Thus, assuming zero bias error, q k+1 remains same as a priori attitude q − k+1 sincẽ b gyro = 0.
By using the previous equations, the attitude error in quaternion is propagated as (17). In this, the direction of ξ k is rotated by a DCM, f DCM q * ω , and subsequently the bias effect is added by qb.
Since both ξ k and qb are small via a small angle assumption, (17) is converted as below.
− k+1 Eq. (18) shows the propagation model for the attitude error in SO (3) form. Note that the approximated equation consists of the gyroscope output and its bias error, which enables to derive the linearized model easily. Taking partial derivatives on attitude error yields Because conceptual attitude error is employed as the navigation state, the linearized model in (19) is arranged in a simpler form than a conventional attitude-based model. Next, the acceleration in the N-frame is depicted by where, Eq. (20) shows the frame rotation of the acceleration by the DCM matrices. Since the error state is employed, two components are considered in the DCM C n b ; one by the estimated attitudê q k and the other by the attitude error ξ k . Here Cb b means the DCM accounting for the estimated attitude perturbation during short period.
To complete linearization for system model, differentiating velocity and position yield the following equations.
In (21), note that the Jacobian term by the attitude error has a skew-symmetric matrix form because of the DCM Cb b . Moreover, by the small angle assumption on ξ k , the DCM C n b can be further approximated with C n b .
The linearized position model is similarly obtained as shown in (22). All the resulting Eqs. (19), (21) ∼ (22) are used together for propagating the state estimate and covariance in the extended Kalman filter.

B. MEASUREMENT MODEL AND LINEARIZATION
This subsection highlights the contribution of the proposed study. In developing measurement model, a formulation is newly established to incorporate pose-dependent magnetic measurement directly in the INS mechanization. Specifically, the magnetic vector itself is modeled by state estimate and the resulting magnetic field residual is used for measurement update in the integration filter. This approach essentially differs from the related work in [21], where the standalone position and attitude solution through burdensome matching method is used for correcting state estimate. Now, given the state variable defined in previous section, the measurement model is represented by the following equation. Note that the model takes a nonlinear functional form, which relates the predicted measurement in the B-frame to the updated position and attitude in the N-frame.
In (23), as discussed before, the DCMs are further separated into terms of attitude and its error. Note that µ i in (23) corresponds with prediction in (9), which matches with the interpolation result via the 3D grid map and a posteriori position estimate.
Eq. (24) and (25) show the linearization result of (23) with respect to the navigation state. The Jacobian by the position can be further derived by differentiating the interpolation equation in Section 2, which is omitted here for brevity.
It can be also observed that the attitude error causes a skew symmetric form involving C b b in (25). Since (23) is independent of the velocity or bias terms, the observation matrix is arranged as (26), which is consequently used for the correction part of the EKF.
Remark: In practical code implementation of EKF, the Jacobian matrix with respect to position in (26) can be derived through the term-by-term partial differential of the interpolation equation. For example, the magnetic field component given as the interpolation function of position can be differentiated term by term with respect to each variable, i.e., 2u 1 0 differentiated from u 2 u 1 with respect to the u axis. This computes ∂m (uvw) /∂u, which yields the 1x1 element of ∂M/∂P and correspondingly, ∂h (x) /∂P. The overall flow chart of the proposed algorithm is shown in Appendix A.

C. INTEGRATED FILTER FORMULATION
Based on the derived system model, the extend Kalman filter is implemented for estimating the integrated navigation results.
Assumed an indirect filter configuration, error states with 15 navigation state elements are introduced with the following definition. x Thus, the INS propagation model is described as where system matrix is given by Equation (28) is the propagation model of the error of the state. It is straightforward to compute the propagation model of error state, where the system matrix F and the input matrix G are already derived. In general, the error covariance can be computed through the system model, although explicit error propagation cannot be acquired owing to its random characteristics. Based on the linearized model, a typical covariance equation is given as where P k ≡ E x k ·x T k . In (29), it is seen an off-line computation can be achieved with the initial covariance P 0 and the input variance Q. Q is determined by the IMU performance characteristics such as random noise ζ and bias stability η in (28).
Measurement update equation is completed by computing the following a posteriori update of Kalman gain and covariance, as shown in (30). where, It is observed that a solid Kalman filter formulation is adapted in (28), as is supported by the linearized observation matrix H in (24). In this formulation, the measurement variance R essentially determines the characteristics of the integration filter. In this study, variance is tuned to be proportional to the strength of magnetic field as appeared in R.
Remark: Physically, tuning R highly depends on the strength of the magnetic field. In view of the signal-to-noise ratio, the measurement variance gets smaller as the magnetic field strength gets larger. Assuming a fixed strength, however, the larger strength implies the closer region around electromagnet. In the region of strong field, due to steep gradient in magnetic field intensity, interpolation error effect becomes more dominant than the signal/noise ratio. In this background, the measurement variance is set to be proportional to the magnetic strength. Note that this is equivalent with the cost function design for pose matching in the previous study [21].
Finally, navigation state is updated by measurement correction; T . Note that Z i implies real magnetic measurements in the B-frame for each carrier, which is obtained by processing magnetometer output as shown in (5).
As observed in (31), a posteriori estimate of navigation state is computed by correcting the measurement residuals. Note that residual between sensor output, z k and,nonlinear prediction via measurement model, h x − k+1 are represented in the B-frame.
The indirect filter yields the result of the error statex in (27). To estimate the navigation state x in (10), we need to linearly add the error state to a priori estimate. Unlike other state variables, however, attitude estimation requires an additional computational step, since attitude error in SO (3) form is assumed in the defined state. To complete the attitude computation, the auxiliary quaternion update follows subsequently as below.
Equation (32) shows that the propagated attitude q − k+1 from (12) is corrected through the attitude error estimate in (31). After correcting attitude, the error k+1 returns to zero for the next propagation step. The following figure shows the overall process of the proposed algorithm.

IV. SYSTEM IMPLEMENTATION
The experimental setup used in this study generally inherits the system concept from previous studies. [21]. While the overall functions of component are similar to the previous work, some improvements are applied to the hardware setup including receiver module and LC resonator part. The receiver sensor module is significantly reduced by substituting electronic devices, passive components and processors as well as optimizing the PCB layout. Through the miniaturized sensor module, application to smaller vehicle is achieved. Figure 6 shows the transmitter system configuration and hardware implementation for generating concurrent AC magnetic fields. In this configuration, 4 exciting carriers are employed with frequencies ranging 40∼70Hz. In detail, the LC resonance circuit consists of a combination of driving transistors, capacitors and inductor coils. For generating excitation current in low frequency ranges, a bundle of capacitors is used in a PCB board while the inductors are installed out of board for mitigating any possible electromagnetic interferences The onboard transistor serves as power converter, which receives an externally generated 40 ∼ 70 Hz digital signal and converts the DC voltage to the corresponding AC voltage. The generated alternating voltage is supplied to the LC circuit and amplified according to the desired resonance frequency to generate a magnetic field. To set the resonance frequency, overall capacitance in the resonance circuit is controlled, while the strength of the AC magnetic field is controlled by the magnitude of the DC voltage applied to the transistor. Consequently, the implemented LC resonance circuit for exciting electromagnets are designed such that it flexibly adjusts field intensity and carrier frequency. Thus, the AC magnetic field coverage can be adaptively tuned on a demanding purpose.  Figure 7 shows the newly manufactured receiver module and onboard package. Compared with the previous sensor module, volume is reduced by 1/4, thus application to various target vehicle can be allowed. The essential magnetic sensor onboard the receiver module is a wide bandwidth, analog magnetometer, the HMC1053 manufactured by Honeywell Inc. Despite its small analog output voltage vulnerable to noise signal, HMC1053 provides fast response time allowing sufficient computational bandwidth.
In our implementation, the noise sensitive magnetic signals are first amplified and filtered via Op-amp, then converted into digital signal via ADC of DSP processor. In the receiver module, additionally, 3D accelerometer and gyroscope are mounted for an inertial integration purpose. Specifically, Murata SCA3100-D07 and InvenSense MPU-3300 are adapted as the onboard accelerometer and the gyroscope, which transmits data as SPI and I2C digital signals, respectively. Along with the converted magnetometer data, inertial sensor measurements are all acquired in the DSP processor, TMS28377D ZWT by Texas Instrument Co., and stored in a micro SDHC card with 1kHz sampling rate. For convenience, storage command transmission and module status check are realized in real time via Bluetooth communication. The 3D map representing the experimental space is modified from the previous research to reflect the characteristics of the new receiver module. This is because the size of new receiver module is reduced significantly. Simultaneously, output sensitivity is changed since the analog signal processing elements are revised in the new PCB. Thus, measurement data for each carrier frequency are compensated by scaling grid map data adaptively. Figure 8 shows an experimental setup of the AC magnetic field environment and operating with the newly manufactured sensor module.

V. EXPERIMENTAL RESULT AND DISCUSSION
The contents of this section are mainly composed of two parts. First, in order to verify the validity of the proposed tightly coupled integration strategy, the pose estimation result is compared with the previous method of map-matching method is presented. Then, the pose estimation performance in a three-dimensional space is analyzed using the newly developed sensor module and the proposed integration method.

A. POSE ESTIMATION COMPARISON
In order to perform the performance analysis according to the algorithm implementation method, a comparative study with the author's previous method was presented. For an objective comparison, the sensor data obtained in the same experiment is used. As mentioned, the comparative algorithm presented in the previous work [21] can be regarded as a loosely coupled integration method since it first computes VOLUME 8, 2020 the optimal position and attitude estimate through minimizing the objective function in the grid map, then combines the pose estimate with the inertial measurement. On the other hand, the proposed method in this paper directly develops magnetic field representation model in terms of the defined filter state, which further constructs filter measurement structure. Therefore, the proposed scheme takes an approach of a tightly coupled integration. In this case, since the magnetic field representation model can be established irrespective of the actual grid map data, the mutual performance comparison is done using the experimental data acquired in the previous paper. First, as already shown in the previous work, Fig. 9 draws the circle trajectory on the two-dimensional plane and then compares each estimation result using the inertial and magnetic field measurements. In this case, black solid line represents the reference trajectory, and blue solid lines and cyan dotted lines represent the proposed tightly coupled integration (denoted by TC, hereafter) and the map matching-based integration (denoted by LC, hereafter), respectively.
The overall error results show that the two estimation algorithms demonstrate similar error characteristics in both horizontal and vertical accuracy. Observing the error graph in the figure, it can be seen the result of TC is slightly improved at the horizontal axis error in the direction of 1 o'clock and 7 o'clock. In the case of the vertical error, both techniques converge to the reference altitude within sub-centimeter level, but due to the inherent error of the measured magnetic field map, a position error of about 2 cm occurs in partial periods. Figure 10 shows the velocity estimation results using the same sensor measurements. As similar with the position estimation result, the horizontal velocity shows no significant difference in error distribution between TC and LC integration. On the other hand, since this experiment employs a horizontal plane trajectory aligned with the reference coordinate system, the vertical axis velocity should be maintained around zero during test periods. When observing the vertical velocity results, it is shown the TC estimation result has better converging characteristics than LC result. This is because the TC structure adopts an extended observation matrix based on multiple magnetic field vectors, which enhances the observability and provides relatively low error sensitivity to the grid map uncertainty. Therefore, more sensitive error correction is expected in estimating state variables such as accelerometer bias, which is supported by case analysis of velocity error. Finally, attitude estimation result is shown in Fig.11. On contrast to position and velocity result, the LC estimation result shows slightly better characteristics in error drift than TC result. Despite an enhanced observability of TC integration, magnetic measurement in relatively remote region is vulnerable to signal fading effect due to low signal-to-noise ratio, thus estimation performance for states with virtually stationary dynamics can be degraded.
Besides, slight offset from reference horizontal plane originates from the mounting offset of sensor platform.
However, in yaw estimation, as large dynamic motion is applied, a negligible estimation difference is observed. In addition, when comparing the average of a 10-second static period at start section and end section, the performance of TC guarantees better repeatability. In addition to accuracy, error distribution can characterize the reliability of integrated navigation system. Figure 12 presents error density graph of each integration method. As discussed previously, magnetic measurement model is primarily based on translational states like position and position increment, accelerometer bias is more sensitive via measurement update while gyroscope bias is less affected. With the enhanced sensitivity to position related states, TC method achieves better error accumulation feature in error range below 0.5cm. While slight improvement is observed in view of estimation accuracy and error distribution, a significant benefit of the proposed TC integration scheme comes from the computational efficiency, which realizes a realtime operation in yielding continuous pose information. Figure 13 compares the computational time of LC and TC. Since the LC technique first solves the map matching problem by searching the grid map iteratively around candidate grids, the required interval for pose estimation is relatively large. On the other hand, TC technique incorporates a direct filter observation model for sequential update whenever magnetic measurement is available, thus greatly reduces computational time. It is remarkable that the total processing time of the TC integration is 35.9 sec, which is shorter than the overall test duration of 40 sec. Consequently, the essential computing advantage of the proposed algorithm facilitates its realtime application onboard a dynamic vehicle.

B. 3D NAVIGATION PERFORMANCE
The quantitative performance of the proposed algorithm including real time capability is presented in the previous section. This section presents the state estimation performance in 3D space by extending the navigation trajectory presented in the previous study. In the 3D trajectory experiment, the newly manufactured miniaturized sensor module mounted on a toy train was used to acquire data. Figure 14 shows the estimated three-dimensional position results for horizontal trajectory in 5 different altitudes. In Fig. 14, the red solid line in the 3 rd plane coincides with the trajectory in the previous section, while other planes are horizontally aligned with only altitude gap by 0.1m. For a detailed analysis, error plots of horizontal and vertical axes are drawn for each plane in Fig. 15 and Fig. 16. In the figure, red solid line represents the reference trajectory, and magenta dot represents the estimation results using the map matching method (Matching), green and blue solid line represents the loosely coupled integration(LC) and the proposed tightly coupled integration(TC), respectively. When observing the figure, the estimation result about the trajectory of the center plane reveals to be most accurate, whereas the error tends to increase as the altitude difference occurs. This is because the magnetic field transmitter is aligned with the third plane, and the accuracy of the magnetic field representation model using the interpolation method gradually decreases as the elevation changes.
The positioning accuracy characteristic in accordance with the interpolation model can be referred to author's previous paper. In addition, in the 1st, 2nd, and 4th planes, it is possible to identify narrow intervals in which the position estimation result is noticeably separated from the reference trajectory. Since the error characteristic does not accumulate and recovers the normal estimation performance after intervals, it is  recognized as the implementation error while constructing the magnetic grid map. Table 1 summarizes the position estimation performance for the three-dimensional trajectory experiment. Root mean square (RMS) error is computed for each floor of Fig. 14. From the results in Table 1, it is observed that the LC and TC have similar error characteristics in the horizontal plane within 0.1 cm. On the other hand, in the case of the vertical axis, the accuracy of the TC is slightly better than about 0.2 ∼ 0.5cm level. In principle, the two estimation schemes demonstrate sub-cm accuracy and show similar estimation performance. It is notable more significant error factor originates from implementation error during sensor map generation as shown. This can be observed in Fig. 15 and 16 as the positioning result via map matching (denoted by MPS) displays amplified variance distributions over some operation periods.
In summary, the proposed algorithm yields similar or slightly better accuracy than the previous map matching method even if considered are the uncertainty of the grid map and the varied filter optimization. Comparatively, the computational time is drastically shortened such that a real-time performance can be achieved. The overall estimation time approximates to 180 sec while the real running time of 5 rotations takes 200 sec in the experiment.

VI. CONCLUSION
This paper proposes a system and algorithm for performing three-dimensional integrated navigation using concurrent AC magnetic fields. In particular, as compared to the author's previous study, the real-time applicability has been greatly improved by implementing a significantly reduced computational time, even with similar estimation accuracy. The main principles for securing the proposed performance improvement are as follows. First, an efficient measurement model between the AC magnetic field vectors and the navigation state was derived by adapting the three-dimensional spatial representation model. Next, the attitude error model was introduced in the integrated filter to secure efficiency in terms of navigation state variable reduction and linearization. A new miniaturized sensor module with faster processing capacity was manufactured to verify the proposed algorithm. Through the experiment in the laboratory environment, it was proved that the accuracy of sub-cm and the real-time computation capability are secured.
In the future, we plan to investigate continuous navigation algorithm for pose estimation through selective integration of AC magnetic fields in the extended test space. For this purpose, it is expected the presented measurement model-based integration suggests a theoretical framework.

APPENDIX
A flow chart that implements the extended Kalman filter with the proposed inertial-magnetic integration is suggested as below. JAEHYU BAE received the B.S. and M.S. degrees in aerospace information engineering from Konkuk University, Seoul, South Korea, in 2017 and 2019, respectively, where he is currently pursuing the Ph.D. degree. His research interests include automatic guidance, navigation and control of aircraft, and unmanned vehicle systems. VOLUME 8, 2020 YOUNG JAE LEE received the B.S. and M.S. degrees in aerospace engineering from Seoul National University, in 1982 and 1985, respectively, and the Ph.D. degree from The University of Texas at Austin, in 1990. He is currently a Professor with the Department of Aerospace Information Engineering, Konkuk University, Seoul, South Korea. His research interests include multiconstellation GNSS system application, integrity monitoring of GNSS, GBAS, RTK, attitude determination, orbit determination of LEO satellites, and integrated navigation-related engineering problems.
SANGKYUNG SUNG (Member, IEEE) received the B.S., M.S., and Ph.D. degrees from Seoul National University, Seoul, South Korea, in 1996Korea, in , 1998Korea, in , and 2003, respectively, all in electrical engineering. From March 1996 to February 2003, he was with the Automatic Control Research Center, Seoul National University. Before joining Konkuk University, Seoul, as an Assistant Professor, in 2007, he was a Senior Researcher with the Telecommunication Research Institute of Samsung Electronics, Suwon. He is currently a Full Professor with the Department of Aerospace Information Engineering, Konkuk University. His research interests include avionic system design and experimentation, inertial and integrated navigation, filter and sensor fusion, urban and seamless navigation, and application to unmanned system. His awards and honors include the Annual Excellent Paper Award (Korean Federation of Science and Technology Societies), the Annual Research Prize (ICROS, KSAS), the Outstanding Conference Paper Award (KIEE), the GSD Innovation Award (National Instrument Company), and the Best Student Paper Award NAECON 2000, the IEEE Aerospace and Electronic Systems Society.