External Disturbances Rejection for Vector Field Sensors in Attitude and Heading Reference Systems

The attitude and heading reference system (AHRS), which consists of tri-axial magnetometer, accelerometer, and gyroscope, has been widely adopted for three-dimensional attitude determination in recent years. It provides an economical means of passive navigation that only relies on gravity and geomagnetic fields. However, despite the advantages of small size, low cost, and low power, the magnetometer and accelerometer are susceptible to external disturbances, such as the magnetic interference from nearby ferromagnetic objects and current-carrying conductors, as well as the motional acceleration of the carrier. To eliminate such disturbances, a vector-based parallel structure is introduced for the attitude filter design, which can avoid the mutual interference between gravity and geomagnetic vectors. Meanwhile, an approach to estimate and compensate the external disturbances in real time for magnetometer and accelerometer is also presented. Compared with existing designs, the proposed filter architecture and external disturbance rejection algorithm can feasibly and effectively cooperate with mainstream data fusion techniques, including complementary filter and Kalman filter. According to experiment results, in the case that large and persistent external disturbances exist, the proposed method can improve the accuracy and robustness of attitude estimation, and it outperforms the existing methods such as switching filter and adaptive filter. Furthermore, through the experiments, the critical role of fading factor in handling the external disturbance is revealed.


Introduction
Attitude and heading reference systems (AHRS) can provide three-dimensional attitude information [1][2][3][4][5][6], and they are widely used in unmanned aerial vehicles (UAV), mobile robots, motion tracking, etc. A typical sensor configuration in AHRS consists of a tri-axial magnetometer, a tri-axial accelerometer, as well as a tri-axial gyroscope. The combination of magnetometer, accelerometer, and gyroscope is also referred to as a MARG sensor [7][8][9], which has the advantages of ultra-low size, cost, and power.
Three-dimensional attitude estimation in AHRS mainly relies on two natural vector fields, namely gravity and geomagnetic fields [1,10]. The former points to the center of the Earth and can be measured by the accelerometer, while the latter points to the magnetic north and can be measured by the magnetometer. Besides the above two vectors, the gyroscope can provide the angular velocity of the carrier that the AHRS is attached to, and it can help to augment dynamic attitude accuracy.

MARG Sensor Error Modeling
As stated above, the MARG sensor in AHRS are used to measure three different vectors, namely the gravity vector g, the geomagnetic vector h, and the angular velocity ω. The measurement model of MARG sensor can be written as (1), which includes various error sources [1, [11][12][13][14][15][16][17].
v acc = C s f , acc C ma,acc C no,acc g + f a + b acc + ε acc v mag = C s f ,mag C ma,mag C no,mag C si (h + b hi ) + b mag + ε mag v gyr = C s f ,gyr C ma,gyr C no,gyr ω + b gyr + δb gyr + ε gyr (1) In Equation (1), the subscripts 'acc', 'mag', and 'gyr' indicate the parameters corresponding to accelerometer, magnetometer, and gyroscope, respectively. The 3 × 1 vector v, b, and ε denote the sensor outputs, biases, and noise terms, respectively. Moreover, the 3 × 3 matrices C s f , C ma , and C no Micromachines 2020, 11, 803 3 of 15 stand for the scale factors, misalignment, and non-orthogonality, respectively. Meanwhile, there are several noticeable issues in Equation (1) for each sensor.
According to Equation (1), the accelerometer is sensitive to both the specific force f a (defined as the non-gravitational force per unit mass) and gravity g. In the case that accelerometer is used to measure the motional acceleration [38], the measurand is f a , while gravity g can be viewed as the disturbance. However, in AHRS, g is needed for attitude estimation, and f a plays the role of external disturbance.
For the magnetometer, the magnetic interferences can be categorized into hard-iron and soft-iron disturbances. The soft-iron disturbance is proportional to the external magnetic field, and it is described by the matrix C si in Equation (1). On the other hand, the hard-iron disturbance usually comes from the permanent magnetism of nearby ferromagnetic materials, and it is described by the 3 × 1 vector b hi in Equation (1). In AHRS, the geomagnetic vector h is used for attitude estimation, and b hi is the source of external disturbance.
For the gyroscope, δb gyr denotes the drift of its bias b gyr , and it should be properly handled in attitude estimation algorithm to avoid accumulative error.
The measurement model in Equation (1) can be rewritten as Equation (2), in which g * , h * , and ω * stand for the measurements of g, h, and ω, respectively. Moreover, the 3 × 3 matrix K and 3 × 1 vector b with corresponding subscripts indicate the deterministic sensor errors. Moreover, δb gyr is simplified to δb, while d a and d m indicate the undetermined external disturbances for the accelerometer and magnetometer, respectively.
The commonly used calibration methods for MARG sensor (such as the ellipsoid fitting method and multi-position method) can determine K and b, but most of them presume that both K and b are time-invariant. In other words, such calibrations are based on linear time-invariant (LTI) error models, and thus they are insufficient to compensate δb, d a , and d m .
In the following discussion, it is presumed that the time-invariant error terms of MARG sensor have already been calibrated appropriately.

MARG Sensor Data Fusion
As stated above, KF and CF are the two major types of MARG sensor fusion techniques. There are two essential parts in KF, namely the state transition (or prediction) and state correction steps. As can be seen in Figure 1a, when KF is used for attitude estimation, in each time step, the 3D attitude Θ is predicted according to ω * , and then the prediction is corrected according to g * and h * . The superscripts '−' and '+' indicate the a priori and a posteriori estimation, respectively.
On the other hand, complementary filters perform data fusion in frequency domain, as shown in Figure 1b,c. In Figure 1b, H(s) stands for a high-pass transfer function that can help to eliminate the gyro bias drift, and TRIAD refers to the tri-axial attitude determination algorithm [43]. Figure 1c shows another widely used implementation of CF that first presented in [18].
Hereafter, it is assumed that the noise term ε and the gyro bias drift δb can be properly handled by the attitude filter (i.e., KF or CF), and thus the only problem left is to handle the external disturbances d a and d m . Hereafter, it is assumed that the noise term and the gyro bias drift can be properly handled by the attitude filter (i.e., KF or CF), and thus the only problem left is to handle the external disturbances and .

Switching Filters
Different criteria can be used to detect the disturbances and/or :

-
The norm of * (or * ): it is the most widely used criterion, due to its high convenience and feasibility.

-
The deviation between and * (or and * ): it is another frequently used criterion, which can be easily acquired from the residual (or innovation) in the state correction step of KF. However, it may be affected by the filter's instability. - The angular velocity: it can be used to detect rotational acceleration (especially the centripetal acceleration) [28], but it is not always reliable since the turning radius is usually unknown.

-
The magnetic dip (i.e., the angle between and the horizontal plane): it can be used to detect [21,27], but with the requisite that the pitch and roll angles are accurately known.
Besides the variety of criteria, there are also different ways to implement hard-switching and soft-switching filters: (a) Tune the filter gain (or the weighting coefficient, cut-off frequency, etc.) directly [8,27]. (b) Enlarge the measurement covariance of * and/or * , to adjust the filter gain indirectly [1, [23][24][25]. This is only fit for KF. (c) Tune the filter according to more complex switching logic, such as the hidden Markov model in [28] and the state machine in [30].

External Disturbances Estimation
Model-based disturbance rejection method tries to extract the external disturbances from raw measurement, rather than discarding it. Therefore, it can probably preserve more information than the threshold-based methods.
To estimate the external disturbance, it can be included in the state variables of attitude filter (specifically the KF) [31,32]. However, it will increase the dimension of KF, as well as the computational burden.
The principle of exterior estimator can be summarized as follows.
(1) Assuming that the a posteriori estimation of and at the (k−1)th time step are already available (denoted as , and , ), the a priori estimations at the kth time step (denoted

Switching Filters
Different criteria can be used to detect the disturbances d a and/or d m : - The norm of g * (or h * ): it is the most widely used criterion, due to its high convenience and feasibility. - The deviation betweenĝ and g * (orĥ and h * ): it is another frequently used criterion, which can be easily acquired from the residual (or innovation) in the state correction step of KF. However, it may be affected by the filter's instability. - The angular velocity: it can be used to detect rotational acceleration (especially the centripetal acceleration) [28], but it is not always reliable since the turning radius is usually unknown. - The magnetic dip (i.e., the angle between h and the horizontal plane): it can be used to detect d m [21,27], but with the requisite that the pitch and roll angles are accurately known.
Besides the variety of criteria, there are also different ways to implement hard-switching and soft-switching filters: (a) Tune the filter gain (or the weighting coefficient, cut-off frequency, etc.) directly [8,27]. (b) Enlarge the measurement covariance of g * and/or h * , to adjust the filter gain indirectly [1, [23][24][25]. This is only fit for KF. (c) Tune the filter according to more complex switching logic, such as the hidden Markov model in [28] and the state machine in [30].

External Disturbances Estimation
Model-based disturbance rejection method tries to extract the external disturbances from raw measurement, rather than discarding it. Therefore, it can probably preserve more information than the threshold-based methods.
To estimate the external disturbance, it can be included in the state variables of attitude filter (specifically the KF) [31,32]. However, it will increase the dimension of KF, as well as the computational burden.
The principle of exterior estimator can be summarized as follows.
(1) Assuming that the a posteriori estimation of d a and d m at the (k − 1)th time step are already available (denoted asd m,k , the measurements g * k and h * k are corrected according to Equation (4).
(3) Once the a posteriori estimationĝ + k andĥ + k are provided by KF, the estimations of d a and d m can be updated according to Equation (5).
The exterior estimator described by Equations (3)- (5) has been proven to be effective in [33][34][35][36][37], but it can only work with KF. Meanwhile, it still has two major problems that need further explanation.
The first problem is the interpretation of Equation (3), which was called the first-order auto-regressive process [36] or Markov chain process [37] that driven by low-pass filtered white noise [33][34][35]. However, in most cases, the external disturbances d a and d m are unpredictable, i.e., there is no sufficient a priori knowledge about their arising and changes over time. Therefore, Equation (3) is by no means rigorous, and it is certainly not driven by white noise (as presumed in [33][34][35][36][37]).
The second problem is the measurement covariance in KF when using the above exterior estimator. Theoretically, if the covariance ofd − a,k is R a , while the covariance of g * k is R g , the covariance of g c k should be R g + R a . Unfortunately, R a remains unknown in most cases due to the lack of a priori knowledge. In [33][34][35][36][37], R a was calculated according to d a itself. For instance, R a can be simply defined as a scalar matrix in Equation (6) [33,37]: Alternatively, R a can be calculated by a windowed smoothing approach, as in Equation (7) [34]: Obviously, Equations (6) and (7) will enlarge the measurement covariance of g c k once the external disturbance arises, but they are not rigorous either.
In the following section, the versatility of the exterior estimation algorithm will be expanded, and a different interpretation of Equation (3) will be given. Furthermore, it will be proved that the above modification of measurement covariance are unnecessary (and can even bring in negative effects) by experiments in Section 4. Figure 2 shows a universal architecture for attitude estimation as well as external disturbances rejection. This architecture can cooperate with all variants of KF and CF, since the estimation and compensation of either d a or d m is performed outside the filter. Moreover, this architecture incorporates several features that can help to simplify the structure and/or enhance the robustness.

Attitude Filter Architecture
First, the 3D attitude information is preserved in g and h, instead of the commonly used quaternion or DCM. Thus, there is no need to calculate the quaternion or DCM according to g and h, or vice versa. Once g and h have been estimated, the heading, pitch, and roll angles (denoted as ψ, θ, and = ϕ, respectively) can be solved by TRIAD algorithm [43]. This feature is known as the vector-based or sensor-based design [41,42].
Micromachines 2020, 11, x 6 of 15 compensation of either or is performed outside the filter. Moreover, this architecture incorporates several features that can help to simplify the structure and/or enhance the robustness. First, the 3D attitude information is preserved in and , instead of the commonly used quaternion or DCM. Thus, there is no need to calculate the quaternion or DCM according to and , or vice versa. Once and have been estimated, the heading, pitch, and roll angles (denoted as , , and , respectively) can be solved by TRIAD algorithm [43]. This feature is known as the vector-based or sensor-based design [41,42].
Secondly, since the filter is vector-based, it is possible to use two sub-filters to estimate and separately. This is called the parallel filter [37], which can avoid the mutual interference between and , and thus help to restrain the impact of external disturbances. Meanwhile, it is also beneficial for reducing the dimension of each sub-filter.
Last but not least, the external disturbances and are estimated outside the attitude filter, i.e., it is an exterior estimator of and , as introduced in Section 2.4. Thus, it is possible to use different algorithms in the filter, including EKF, UKF, CKF, and CF. In each time step, the proposed attitude estimation and external disturbances rejection algorithm can be summarized as follows:


Step 1: Predict the disturbances , and , . This step is the same as (3).  Step 2: Subtract , and , from the corresponding measurements (i.e., * and * ) to get the corrected vectors (denoted as and ). This step is the same as (4). Secondly, since the filter is vector-based, it is possible to use two sub-filters to estimate g and h separately. This is called the parallel filter [37], which can avoid the mutual interference between g and h, and thus help to restrain the impact of external disturbances. Meanwhile, it is also beneficial for reducing the dimension of each sub-filter.
Last but not least, the external disturbances d a and d m are estimated outside the attitude filter, i.e., it is an exterior estimator of d a and d m , as introduced in Section 2.4. Thus, it is possible to use different algorithms in the filter, including EKF, UKF, CKF, and CF. In each time step, the proposed attitude estimation and external disturbances rejection algorithm can be summarized as follows:

External Disturbances Estimation
• Step 1: Predict the disturbancesd − a,k andd − m,k . This step is the same as (3).

•
Step 2: Subtractd − a,k andd − m,k from the corresponding measurements (i.e., g * k and h * k ) to get the corrected vectors (denoted as g c k and h c k ). This step is the same as (4).

•
Step 3: Get the estimated vectorsĝ k andĥ k according to g c k , h c k , and ω * k through a certain algorithm (KF or CF).

Specific Implementations
The above algorithm can cooperate with KF and CF, as presented below.

The Proposed Method with Kalman Filter (Abbreviated as PM-KF Hereafter)
• PM-KF Step 1: Get a priori estimation according to Equation (8), in which the state transition matrix is calculated as satisfies ω * k × u = ω * k × u for any 3 × 1 vector u, and ∆t is the sampling period.
• PM-KF Step 3: Get a posteriori estimation according to Equations (10)- (12). The measurement covariance matrices are defined as R g = σ g I 3×3 and R h = σ h I 3×3 , with σ g and σ h denoting the noise of accelerometer and magnetometer, respectively.
Then k = k + 1 and go to PM-KF Step 1.

The Proposed Method with Complementary Filter (Abbreviated as PM-CF Hereafter)
• PM-CF Step 1: Correct the angular velocity according to Equation (14).
• PM-CF Step 2: Update gravity and geomagnetic vectors using to the corrected angular velocity, as described by Equations (15)- (17).
Then k = k + 1 and go to PM-CF Step 1.

Remarks
As can be seen in the above implementations, the proposed parallel vector-based architecture can greatly simplify the design of attitude filter, especially the KF. There is no need to use UKF, CKF, or other nonlinear version of KF, since both the process and measurement models are concise and straightforward. Moreover, the proposed external disturbance rejection method seems to be a direct extension of Equations (3)-(5). However, it can be interpreted in a different way as follows.
As pointed out in the above section, the external disturbances d a and d m are unpredictable in most cases. As a matter of fact, the essence of Equation (3) is the assumption that either d a or d m will not have drastic change during a single sampling period, and thus the last a posteriori estimation can be used to correct the sensor measurement at the present time step. Obviously, this approach is inexact, and the actual role of the coefficients c a or c m is to make the estimation error tend to decrease. In other words, c a and c m are actually the fading factors that help to avoid divergence. Nonetheless, if c a (or c m ) is too close to zero, it will inevitably weaken the compensation of d a (or d m ). The effects of these two fading factors will be evaluated in the following section.
To find a general solution for not only the KF, it is unnecessary to modify the covariance of g c k and/or h c k . In fact, the proposed algorithm makes no change to the attitude filter itself, and thus it is an actual exterior estimator of external disturbances. Moreover, it can be seen in the following experiments that the modification of the measurement covariance is not always helpful to improve the performance.

Basic Settings
The proposed algorithm is evaluated on an AHRS module, which is based on a monolithic MARG sensor MPU9250 and working at the sampling rate of 20 Hz. Meanwhile, a single-axis rate table is used to provide heading and angular rate references with 0.0001 • resolution. The AHRS module and rate table, as well as the hardware installation, are shown in Figure 3.
The magnetometer and accelerometer in MPU9250 are calibrated using a calibration scheme based on dual inner products, which can be viewed as an improved ellipsoid fitting method and was detailed in [44]. Moreover, the gyroscope in MPU9250 is calibrated using a cross product-based algorithm, which was presented in [45]. Micromachines 2020, 11, x 9 of 15 The magnetometer and accelerometer in MPU9250 are calibrated using a calibration scheme based on dual inner products, which can be viewed as an improved ellipsoid fitting method and was detailed in [44]. Moreover, the gyroscope in MPU9250 is calibrated using a cross product-based algorithm, which was presented in [45].
After the above calibration, raw data are then acquired according to the flowchart in Figure 4. It can be seen that the AHRS experiences high-speed stop-and-go rotations, which can generate considerable centripetal acceleration. Meanwhile, it also suffers artificial magnetic interference that caused by ferromagnetic objects. Three different algorithms are used for attitude estimation, including PM-KF and PM-CK that introduced in Section 3.3, as well as the parallel Kalman filter in [37] (abbreviated as PA-KF). The After the above calibration, raw data are then acquired according to the flowchart in Figure 4. It can be seen that the AHRS experiences high-speed stop-and-go rotations, which can generate considerable centripetal acceleration. Meanwhile, it also suffers artificial magnetic interference that caused by ferromagnetic objects. The magnetometer and accelerometer in MPU9250 are calibrated using a calibration scheme based on dual inner products, which can be viewed as an improved ellipsoid fitting method and was detailed in [44]. Moreover, the gyroscope in MPU9250 is calibrated using a cross product-based algorithm, which was presented in [45].
After the above calibration, raw data are then acquired according to the flowchart in Figure 4. It can be seen that the AHRS experiences high-speed stop-and-go rotations, which can generate considerable centripetal acceleration. Meanwhile, it also suffers artificial magnetic interference that caused by ferromagnetic objects. Three different algorithms are used for attitude estimation, including PM-KF and PM-CK that introduced in Section 3.3, as well as the parallel Kalman filter in [37] (abbreviated as PA-KF). The Three different algorithms are used for attitude estimation, including PM-KF and PM-CK that introduced in Section 3.3, as well as the parallel Kalman filter in [37] (abbreviated as PA-KF). The main difference between PM-KF and PA-KF is that the latter has adaptive measurement covariance, as described by Equation (6).
In both PM-KF and PA-KF, the sensor noise covariances are σ h = 0.5 µT, σ g = 0.1 m/s 2 , and σ ω = 0.005 rad/s, respectively. On the other hand, the coefficients of PM-CF are set to K g = K h = 0.3.

Experiment Results
Since each of the above algorithms outputs the estimations of g and h, their performance is evaluated in terms of the directional errors of g and h, i.e., the angle between the estimated vector and the true vector. At the kth time step, if the estimated vector isĝ k (orĥ k ), while the true vector is g r,k (or h r,k ), the directional error can be calculated as δ g,k = cos −1 ĝ k ·g r,k / ĝ k g r,k (or δ h,k = cos −1 [(ĥ k ·h r,k )/( ĥ k h r,k )]). Moreover, both fading factors c a and c m are increased from 0 to 2 with the increment of 0.01, so as to evaluate their impacts. Figure 5 shows the maximum error and root mean square error (RMSE) of each algorithm versus the fading factors. main difference between PM-KF and PA-KF is that the latter has adaptive measurement covariance, as described by Equation (6).
In both PM-KF and PA-KF, the sensor noise covariances are 0.5 μT, 0.1 m s ⁄ , and 0.005 rad s ⁄ , respectively. On the other hand, the coefficients of PM-CF are set to 0.3.

Experiment Results
Since each of the above algorithms outputs the estimations of and , their performance is evaluated in terms of the directional errors of and , i.e., the angle between the estimated vector and the true vector. ). Moreover, both fading factors and are increased from 0 to 2 with the increment of 0.01, so as to evaluate their impacts. Figure 5 shows the maximum error and root mean square error (RMSE) of each algorithm versus the fading factors. As can be seen in Figure 5, the fading factors and have significant impacts on the performance of each algorithm. It is clear that all the algorithms will diverge once or is greater than 1. Meanwhile, in Figure 5, the optimal fading factors for each algorithm (i.e., the fading factors that lead to the best performance) can be found, as listed in Table 1. It is noteworthy that the optimal fading factors of PA-KF are quite small, especially the optimal value of is only around 0.05 (but coincides with [33] and [37]). As mentioned in Section 3.4, such small fading factors will weaken the compensation of external disturbances.
To further demonstrate the performance of each algorithm with its optimum fading factors, the estimated vectors and , the estimated disturbances and , as well as the corresponding 3D attitude outputs (heading , pitch , and roll ) are all plotted in Figures 6 and 7. As can be seen in Figure 5, the fading factors c a and c m have significant impacts on the performance of each algorithm. It is clear that all the algorithms will diverge once c a or c m is greater than 1. Meanwhile, in Figure 5, the optimal fading factors for each algorithm (i.e., the fading factors that lead to the best performance) can be found, as listed in Table 1. It is noteworthy that the optimal fading factors of PA-KF are quite small, especially the optimal value of c a is only around 0.05 (but coincides with [33] and [37]). As mentioned in Section 3.4, such small fading factors will weaken the compensation of external disturbances. To further demonstrate the performance of each algorithm with its optimum fading factors, the estimated vectors g and h, the estimated disturbances d a and d m , as well as the corresponding 3D attitude outputs (heading ψ, pitch θ, and roll ϕ) are all plotted in Figures 6 and 7.
In Figure 6a,b, the centripetal acceleration caused by rapid rotations can be spotted, and it results in the fluctuation and divergence in the output of PA-KF in Figure 6d,e. On the contrary, the outputs of PM-KF and PM-CF remains stable when the centripetal acceleration exists. This is consistent to the results in Table 1, i.e., both PM-KF and PM-CF have better results than PA-KF in compensating the large and lasting external acceleration.
Nevertheless, it can also be noticed in Table 1 that PM-KF and PM-CF do not outperform PA-KF when compensating magnetic interference. More details of the magnetic interference and its compensation by different algorithms are shown in Figure 7.

Discussion
The above experiment results have proven the feasibility and effectiveness of the proposed method, including the vector-based parallel architecture and the external disturbance rejection algorithm. Still, there are some noteworthy issues.
First, the proposed method (including PM-KF and PM-CF) outperforms PA-KF when dealing with the centripetal acceleration, but shows no superiority in handling the magnetic interference. As shown in Figure 6, the raw data from accelerometer contain large and long-term centripetal acceleration, and the proposed method is more adaptive to such hostile conditions. On the contrary, PA-KF works better against momentary and moderate disturbances, as shown in Figure 7.
The main reason for the above phenomenon lies in the data fusion process of MARG sensor. Once the external disturbance arises, the measurement of accelerometer and/or magnetometer becomes unreliable, and thus the attitude estimation algorithm can only rely on the integration of angular velocity. However, inevitable bias drift of gyroscope results in a dilemma that the attitude filter will either diverge due to accumulated error, or still use the distorted measurement of accelerometer and/or magnetometer. In a word, any solution for disturbance rejection is essentially a trade-off between the gyro drift and the external disturbances.
Compared to the switching filters that completely discard the unreliable measurements, both PA-KF and the proposed method partially retain the measurements of accelerometer and magnetometer even when they contain external disturbances. However, PA-KF assigns more weight to the gyroscope by enlarging the measurement covariance of accelerometer and/or magnetometer. Consequently, PA-KF will be more significantly affected by the gyro drift when the disturbance is large and lasting.
Another key issue is the fading factor c a and c m , which can greatly impact the compensating of external disturbance. As discussed in Section 3.4, a fading factor between 0 and 1 can help to decrease the estimation error, and it is proven by the fact that the optimal fading factors of all the algorithms are less than 1 in Table 1. However, if the fading factor is too small, it will evidently weaken the compensation of external disturbance, and it is also proven by the results listed in Table 1.
Compared to PM-KF and PM-CF, the optimal fading factors of PA-KF are much closer to 0. Since PA-KF is more significantly affected by the gyro drift due to its adjustment to the measurement covariance, it needs a smaller fading factor to alleviate this problem. Unfortunately, such a small fading factor will definitely weaken the compensation effect, and it results in the poor performance of PA-KF when handling large and lasting centripetal acceleration.
Finally, it is worth mentioning again that the main advantage of the proposed method is not only its better performance against strong and persistent disturbances, but also its flexibility to cooperate with commonly used sensor fusion algorithms (including but not limited to KF and CF).

Conclusions
In this paper, a versatile external disturbance rejection approach is presented, along with a vector-based parallel architecture for 3D attitude estimation. It is proven by experiments that the proposed filter architecture and disturbance rejection approach can work well with KF and CF in the presence of external acceleration and magnetic interference.
The proposed method provides a feasible and generally applicable solution for MARG-based 3D attitude estimation. Nonetheless, the problem of external disturbance rejection for accelerometer and magnetometer is far from completely solved. As stated in the above section, all solutions for disturbance rejection in MARG-based AHRS, either the existing approaches or the proposed method, are no more than trade-offs between the gyro drift and the external disturbance. Since the significant role of the fading factor is revealed, a possible direction for future exploring is the adaptive adjustment of such fading factors.