Reference Frame Unification of IMU-Based Joint Angle Estimation: The Experimental Investigation and a Novel Method

Inertial measurement unit (IMU)-based joint angle estimation is an increasingly mature technique that has a broad range of applications in clinics, biomechanics and robotics. However, the deviations of different IMUs’ reference frames, referring to IMUs’ individual orientations estimating errors, is still a challenge for improving the angle estimation accuracy due to conceptual confusion, relatively simple metrics and the lack of systematical investigation. In this paper, we clarify the determination of reference frame unification, experimentally study the time-varying characteristics of reference frames’ deviations and accordingly propose a novel method with a comprehensive metric to unify reference frames. To be specific, we firstly define the reference frame unification (RFU) and distinguish it with drift correction that has always been confused with the term RFU. Secondly, we design a mechanical gimbal-based experiment to study the deviations, where sensor-to-body alignment and rotation-caused differences of orientations are excluded. Thirdly, based on the findings of the experiment, we propose a novel method to utilize the consistency of the joint axis under the hinge-joint constraint, gravity acceleration and local magnetic field to comprehensively unify reference frames, which meets the nonlinear time-varying characteristics of the deviations. The results on ten human subjects reveal the feasibility of our proposed method and the improvement from previous methods. This work contributes to a relatively new perspective of considering and improving the accuracy of IMU-based joint angle estimation.


Introduction
Estimating angles of human limb joints plays a fundamental role in obtaining human kinestate and assessing human ability, which is involved in various fields ranging from rehabilitation to robotics, from industry to household. Clinicians observe the joint rotation conditions of lower-limb joints to diagnose the ongoing progression of gait-related diseases, such as muscle dystrophy, joint injuries and stroke, etc. [1][2][3][4][5], based on which further treatment like surgical intervention, recovery or rehabilitative therapy [6,7] will be applied. As for robotics, joint angle estimation contributes to the robots' awareness of human kinestate and thus results in a smart and accurate assistive and/or rehabilitative action during physical human-robot interaction scenarios [8,9]. Moreover, the angle estimationinvolved techniques enable the assessment of human ability. For example, studies predicted risk-related information based on the precepted joint angle information, like [10] for fall prevention, ref. [11] for working safety assessment. Other studies utilized joint angles to assess motor-related information, like [12] for metabolic cost estimation, ref. [13] for mobility assessment for sports. Given that joint angle estimation functions importantly in Pioneer studies like [23,28] utilized the consistency of the angular rates of one calibration motion to measure and compensate for the deviation. Following studies like [29,30] leveraged the joint axis of the hinge-joint constraint [27] to online detect the deviations and calculate the compensation. Some other studies [13] employed the consistency of the gravity acceleration. Secondly, the characteristics of the reference frames' deviations should be systematically investigated so as to propose a suitable RFU method. The pioneer study [28] unified reference frames in a static manner. The following study [23] proposed a linear interpolated method that improve the RFU efficiency. Latter studies started to employ a point-wise compensation for RFU [13,29,30]. Although presenting increasing improvement, the compensation methods for RFU are proposed empirically without the instruction of reference frame deviation characteristics.
In this study, we systematically define the RFU problem and carefully distinguish it from the drift correction problem. Then, we analyze the deviations of reference frames caused by the different characteristics of IMUs, in the manner of excluding the confounding factors. Based on the analysis, we propose comprehensive metrics and further utilize it to form a novel RFU method that meets the characteristics of deviations. Both specially devised three-dimensional gimbal and human movement experiments are employed to separately estimate the performance. Performance comparison with previously used methods demonstrates the feasibility and improvement of our proposed method.

Related Work
Before presenting our method, it is important to firstly clarify some issues like what is the RFU problem, what is expected to know about the deviation of reference frames and how we get inspirations from literatures.
As stated above, reference frame unification is to let the sensor-fixed frames of IMUs depicted in the same reference frame, with the aim of decreasing the errors caused by deviations of reference frames and thus improving the accuracy of angle estimation. It should be noted that although most literatures confuse the "drift correction" [29,30,32,33] and the "reference frame unification", RFU is significantly different from drift correction (DC). In our paper, we will systematically define the RFU problem, distinguish it with the DC problem and mathematically formulate them.
Previously used methods of calculating the rotational compensation of RFU are based on different hypotheses of reference frame deviations' characteristics. Pioneer works like [23,28] assumed the time-varying characteristics of the deviations can be approximated statically or linearly. Following works utilized the consistency of the joint axes estimated by the hinge-joint constraint to calculate point-wise rotational compensation [29]. Although the performance is limited by the leveraged single metrics, the point-wise compensation theoretically provides a better approximation to the time-varying deviations. Other studies assumed the deviations can be statically compensated through fusing the corrections calculated by different metrics [31]. In our study, we experimentally study the time-varying characteristics of the deviations, and utilize the conclusions to propose a novel point-wise and time-varying compensation for RFU.
The metrics used to measure the deviations also relate to the calculation of the rotational compensation. Studies [23,28] calculated the rotational compensation based on the consistency of angular rates of IMUs mounted on thigh and shank during the hip abduction/adduction movement. The angular rate-based metric may be sensitive to the movement of other lower-limb joints. For example, if the knee rotated during the hip abduction/adduction movement, the consistency of the angular rates would be broken and would significantly degrade the RFU performance. Moreover, the implementation of calibration motions might impede the angle estimation technique's broader usage. Ref. [29] proposed to measure the deviations by the joint axis of the hinge-joint constraint. Following the hinge-joint constraint, lower-limb joints were simplified as a 1-DoF joint, and the coordinates of the joint axis in sensor-fixed frames can be estimated. Through calculating the rotational deviation of the joint axis' coordinates in different reference frames, the devi-ations of reference frames can be measured and then compensated. The joint axis-based metric might be limited by the hinge-joint constraint that simplifies the biological joints as 1-DoF joints [27]. Although the study [29] has proposed to screen the hinge-rotating case and interpolate the rotational compensation during non-hinge-rotating cases, as argued above, the linear interpolation still remains limitations. Fasel et al. [13] proposed to measure the deviations using the consistency of the gravity acceleration. That is, the accelerations measured by IMUs should process the same coordinates in the same reference frame. This metric, although promising, is limited by the corruption of movement-caused accelerations. In contrast to the single metric, Seel et al. [31] fused the rotational compensation calculated by the consistency of the local magnetic field and the gravity. This more comprehensive metric provided a better solution to decreasing the flaws of each single metrics, but was still impeded by the static fusion coefficient that needs manually tuning and the not-comprehensive-enough metrics. Inspired by the existing metrics and the lessons of their flaws, we propose a novel metric that comprehensively utilizes all the above-mentioned metrics, and thus contributes to a more accurate RFU.

Method and Materials
In this section, we firstly mathematically formulate the problem and secondly describe the experimental procedure and how the experiments are designed according to the paradigms of other works. Then, the experimental investigation is given in order to analyze the characteristics of reference frames' deviations and to instruct the design of our RFU method. Finally, our RFU method is proposed by following the investigated characteristics and principles of reference frames' deviations.

Problem Statement
The RFU and DC problems can be mathematically formulated as follows. The measurements of the ith IMU can be described asω i = ω i + δ ω i ,ã i = a move i + g + δ a i , m i = m i + δ m i , where g denotes gravity, δ denote noise and corruptions that differ across IMUs,ω i ,ã i ,m i denote measurements and ω i , a i , m i denote real values of angular rate, movement-caused acceleration and local magnetic field. The absolute orientations estimated by the corrupted measurements are denoted by the rotation matrices R g 1 s 1 , R g 2 s 2 , which depict sensor-fixed coordinate frames' orientations with respect to different reference frames. It should be noted that in the following, we treat the quaternion q and the rotation matrix R as equal and without stating the substitution between the two again. The RFU issue can be formulated as where q corr denotes the rotational compensation of RFU. The DC issue can be formulated as: where q dri f t denotes the rotational compensation of DC, r s i and r g i denote reference vectors in the sensor-fixed frame and the Earth frame, respectively. RFU can be significantly distinguished from DC from two aspects. On one hand, the result of RFU and DC refers to different coordinate frames. As shown in Figures 1 and 2, RFU is to make the vectors in sensor-fixed frames depicted in the same reference frame, regardless of what the reference frame is. DC is to correct the drift of absolute orientation estimation mainly caused by magnetic distortion, expecting to make the reference frame the same as the Earth frame.
On the other hand, RFU and DC differ by the calculation of the rotational compensation. Due to the magnetic measurements' dominant influence of determining the azimuth angle of orientation, the operation of DC is to mainly calculate the rotational compensation of the horizontal plane. In contrast, RFU is just to rotate one reference frame to correspond to another one, without consideration of the rotational axis.

Experimental Procedure
In order to study the characteristics of different reference frames' deviations, other factors during angle estimation, like sensor-to-body alignment, relative rotation between two segments, need to be decreased as much as possible. A gimbal with flat surfaces and known rotational axes was designed to mimic the biological lower-limb joints, mainly following the paradigm of the design in [23]. As shown in Figure 3, the Hall sensors attached to each axis by couplings were used to measure angles of each rotational axis, the signals of which were sampled at 500 Hz. Four IMUs (74Hz, Delsys Trigno, IM type & Avanti type) were attached to each segment of the gimbal. Only IMU2 and IMU3, which were attached to the segments beside a 1-DoF joint (the joint axis j 2 ), were used in this study. Z axis of IMU2 and X axis of IMU3 were placed in the direction of the main axis j 2 . The null position of the axis j 2 is to let the X-Z plane of IMU2 and Y-Z plane of IMU3 stay in the same plane, which is ensured by the devised mechanical structure of the gimbal. Due to the flat mounting surface, IMUs can be mounted with known orientation relative to segment-fixed frames by the means of manually aligning the sides of IMUs and gimbal following the manners of [23,34]. In this way, errors of the sensor-to-segment alignment can be maximally excluded from the calculation. The influence of relative rotation between two segments, i.e., the mounting surfaces besides the joint axis j 2 , will be excluded by calculation stated in latter sections.  In order to evaluate the superiority of our proposed RFU method, level walking experiments on human subjects were performed. Ten healthy subjects without previous neuropathological history (7 males and 3 females, age range: 20-30 years, height range: 155-184 cm, weight range: 50-90 kg) were recruited and asked to perform five trials of 3-min level walking. Rest periods were allowed between trials to avoid fatigue. The experiment protocol was approved by the local ethical committee and all participants had been informed of the content and their right to withdraw from the study at any time, without giving an explanation.
Before and after the five trials, calibration postures proposed in [25] were asked to be performed by each subject in order to calculate the 3-dimensional joint angles of each joint. Besides, the hip abduction/adduction movement, following the paradigm of [23,28], was also performed in order to calculate the static and linearly interpolated RFU based on the angular-rate metrics. As shown in Figure 4, multiple IMUs were attached to the trunk, right thigh, right shank and right foot. The redundant IMUs were attached to study the influence of IMUs' placement, which will feature in our future study. In this work, we randomly select one IMU per segment without loss of generality. Given that we employ magnetometer measurements in IMUs' orientation estimating, ferromagnetic metal plates were placed in the proximity of IMUs in order to induce local magnetic distortion, which followed the paradigm of [30]. Sixteen retro-reflective markers were attached to subjects' pelvis and lower limbs following the principles of [14,15], whose 3-D locations were recorded (100 Hz) using an 8-camera video system (Vicon, Oxford, UK). The joint angles were calculated by the pose estimation and inverse kinematics model embedded in the software Visual 3D. The signals from nine-axis IMUs and the video system were synchronized by the trigger and time stamps.

Data Preprocessing
In order to evaluate the deviations' characteristics of different IMUs, primary factors like different angular rates and IMUs' individual characteristics should be included, while confounding factors like sensor-to-body alignment, rotation-caused IMUs' orientation difference should be excluded. To be specific, confounding factors can be excluded by the careful alignment of IMUs and the gimbal's mounting surface and some calculations. As shown in Figure 5, on the null position, the Z axis of the sensor-fixed frame of IMU2, the Y axis of the sensor-fixed frame of IMU3 and the axis j 2 are in the same direction, and the X-Z plane of IMU2 corresponds to the Y-Z plane of IMU3. The determination of the sensor-fixed frames of IMU2 and IMU3 is adopted from the user guide of Delsys.
Under this circumstance, the rotation matrix between the sensor-fixed frames of IMU2 and IMU3 can be given as where T I MU3 I MU2 denotes the rotation that transform the orientation of the sensor-fixed frame of IMU2 to that of IMU3. If the axis j 2 is manually rotated away from the null position, as shown in Figure 5, the consequent rotation can be depicted by the matrix where R denotes the rotation matrix, θ denotes the rotational angle around the axis j 2 , measured by the Hall sensor. In so doing, the orientation of the sensor-fixed frame of IMU2 can be transformed into that of IMU3 by the cascade multiplication of T I MU3 I MU2 and R. Thus, confounding factors are excluded. As for the primary factors, the manually activated rotation around j 2 contributes to different angular rates of the two IMUs thus affects the absolute orientation estimating of both IMUs, while employing two different IMUs also contributes to the different IMUs' individual characteristics.  We next implement the absolute orientation estimating methods adopted from [22] in order to estimate the absolute orientations of IMU2 and IMU3. The estimates are denoted by rotation matrices, i.e., R g 2 s 2 for the absolute orientations of IMU2 and R g 3 s 3 for absolute orientations of IMU3. The estimation accuracy denoted by the root mean square error (RMSE) during dynamic movement was reported to be 1.63 deg in [22]. Following the calculations of excluding the confounding factors, the transformed orientation of IMU2 is calculated as follow.R g 2

Analyzing Procedure
An analysis is performed on the orientation deviations of the two IMUs. Firstly, XYZ Euler angle decomposition is performed to transform the estimated rotation matrices,R g 2 s 2 andR g 2 s 3 , into Euler angles, which are denoted by α i , β i , γ i for X, Y, Z axes, i = 2, 3 for IMU2 and IMU3 respectively. Secondly, differences of each Euler angles, denoted by ∆α, ∆β, ∆γ, are calculated for further analysis. Thirdly, overall deviations and time-evolving characteristics are depicted and analyzed by RMSE and statistical analysis. Particularly, the differences obtained in each trial are divided into three subsets, each lasting 40 s. Then the RMSEs of each subset are averaged over all the five trials. One-way ANOVA is employed to assess the results.

Analyzing Results
In the following, we report our primary findings of the deviations and our analyzing results. Then, we discuss the findings related to the deviations, expecting to indicate the design of our method.
As shown in Figure 6, the differences of the estimated absolute orientations of IMU2 and IMU3 are depicted by the differences of XYZ-decomposed Euler angles. The ranges of ∆α, ∆β, ∆γ are 2.071 deg, 6.836 deg and 7.082 deg, respectively. The relative RMSEs are presented in Figure 7. As shown in Figure 7a, the RMSEs for the roll angle of the three subsets averaged over all the trials are presented. There is a significant difference between the RMSEs of the first two subsets. As shown in Figure 7b, the RMSEs for the pitch angle of the three subsets averaged over all the trials are presented. There is a significant difference between the RMSEs of the first two subsets and between those of the last two subsets. As shown in Figure 7c, the RMSEs for the yaw angle of the three subsets averaged over all the trials are presented. There is a significant difference between the last two subsets in term of RMSE.

Findings
Our main findings basically support our statement in the related work session. That is, (1) the deviations non-linearly evolve over time; (2) the deviations between reference frames are significantly different from the drifts caused by magnetic distortion; (3) the deviations significantly differ from IMUs' estimation drift, especially from heading drift.
Firstly, the significant differences existed among RMSEs of each subset indicate the time-varying characteristics of the deviations. Thus, the static method proposed in [28] cannot provide an acceptable compensation for the deviations, i.e., an accurate unification of reference frames, which meet the conclusion of [23].
Secondly, the time-dependent varying of deviations is non-linear. As shown in curves of Figure 6, the deviations do not vary in a linear or pseudo-linear way. Rather, a combined partially piecewise linear and non-linear manner is presented. This phenomenon suggests that the deviations can be compensated in a piecewise linear way.
Thirdly, heading drift is obviously not the only origin of the deviations. As stated in the related work, most of the related literatures confuse the "drift correction" with the "reference frame unification". That is, the deviations of reference frames are attributed to the heading drift caused by the local magnetic distortion. As shown in Figure 7a,b, the deviations on the pitch and roll angle present a significant amount compared with that of the deviations on the yaw angle. Other than the heading drift that primarily affects the accuracy of yaw, the deviations of the two reference frames, even when confounding factors are excluded, still include estimates' differences on pitch and roll resulted from other factors, such as the movement-caused corruptions of accelerations and measurement noise. Moreover, the RMSEs and range of deviations reported in our analyzing results are relatively larger than those reported in previous IMU orientation estimating studies [21,22]. This also gives side proof of the distinction of "estimation drift" and "reference frame deviation". The larger results of our study result from the combined usage of estimation drift and individual characteristics of IMUs' measurements.
In summary, the findings indicate that the RFU method should not only compensate for the estimates' heading drift but also take other factors into consideration. Moreover, pointwise RFU should be performed to meet the non-linearly time-varying characteristics of the deviations, if possible, with multiple metrics to give a comprehensive measurement.

Reference Frame Unification with Comprehensive Metrics
As stated in the problem statement and Figure 8, the RFU issue is to find and calculate a correction quaternion q corr that rotates [g 2 ] to [g 1 ], thus both absolute orientations are depicted in the same reference frame. Indicated by the findings and literatures, in this section, we propose a novel RFU method that can provide a pointwise compensation with comprehensive metrics. Our method is distinguished from other RFU methods by (1) comprehensively using magnetic fields, accelerations and the coordinates of the main axis as metrics and (2)   In order to calculate q corr , some measurements that are common among reference frames and sensor-fixed coordinate frames should be employed as an intermediate reference to solve the RFU issue. Inspired by literatures, the gravity acceleration, the local magnetic field and the joint axis under the hinge-joint constraint can be leveraged as the intermediate reference. Without considering the movement-caused accelerations, the acceleration vectors, a 1 and a 2 , should be identical in the common reference frame, which contributes to the construction of the correction quaternion, given by where W acc is the rotation axis, θ acc is the rotation angle and q acc is the correction quaternion calculated from gravity. Similarly, the identical magnetic field vectors m 1 and m 2 without magnetic distortion could contribute to another correction quaternion, given by However, due to the corruption of acceleration and magnetometer readings, neither of the correction quaternions calculated above is accurate enough to represent the rotational relationship between [g 1 ] and [g 2 ]. Herein, a weighted sum of these two correction quaternions is used to make a fusion. q corr = k mag · q mag + k acc · q acc (8) In [31], Seel et al. proposed to fuse the two correction quaternions by manually tuning the coefficients. The manually tuned fusion coefficients, as argued in our previous work [22], cannot well fit the time-varying trends of the corruptions of movement-caused accelerations and magnetic distortion. Thus, dynamic tuning should be included. Inspired by the works [29,30] that solely used the joint axes estimated by the hinge-joint constraint, we employ the same paradigm here to form a more comprehensive metric. To be specific, j 1 and j 2 , denoted as the coordinates of the joint axis of a hinge joint depicted in the two sensorfixed coordinate frames, are estimated by the method proposed in [27]. After estimating the absolute orientations, the deviation of [j 1 ] g 1 and [j 2 ] g 1 is introduced here to evaluate to what extent coordinates in [g 2 ] is rotated by q corr into those in [g 1 ], where [j 1 ] g 1 and [j 2 ] g 1 are j 2D 1 and j 2D 2 described in the reference frame [g 1 ],given by where ⊗ denotes the multiplication of quaternions. Regarding f i (k mag , k acc ) = [j 1 ] g 1 − [j 2 ] g 1 as the cost function, fusion coeficients, k mag and k acc , can be estimated through the secant version of L-M method. Hence, the calibration of reference frames synthesizes the information from local magnetic field, gravity and the main axis of the 3-DOF joint, which results in a comprehensive metric. For the sake of real-time calculation, as shown in Figure 9, two sliding windows are separated by an interval (I n ) within which the fusion coefficients estimated in the last sliding window are used to unify reference frames in the interval. At the end of an interval, the last N measurements could be used to compute the cost function and execute iterations to update fusion coefficients.

Data Analysis
Data recorded in the experiment on human subjects are analyzed as follow. Firstly, accuracy is depicted by RMSE between IMU-calculated joint angles and the angles measured by the optical motion capture system, for 3-dimensional angles of hip, knee and ankle. RMSEs of the five trials are averaged, from which the means and standard deviations (SD) are derived in order to evaluate the accuracy and precision. In addition, one-way ANOVA is applied to the performance of different RFU methods. Secondly, in order to evaluate the influence of movements on the correction quaternion, the repeatability test is performed by the dispersion of the five correction quaternions, following the paradigm of [13]. Particularly, the dispersion χ of the five correction quaternions q S,T is calculated around their mean q S for all subjects.

Results
Performance is analyzed following the paradigm in the data analysis section. For the comparison purpose, we implement four previously used RFU methods, which are the static [28] and linearly interpolated [23] methods supported by the hip abduction/adduction before and after measurement, the joint axis-based pointwise RFU method [29] and gravity acceleration-based method [13]. All the data are processed by MATLAB 2015 B.
As shown in Figure 10, all the estimates track the "standard" curves well with similar trends. As shown in Figure 11, our proposed RFU method significantly outperforms other RFU methods (α < 0.05), which are depicted by the means and standard deviations of RMSEs. The RMSEs of each RFU methods averaged over trials and subjects are shown in Table 1. Figure 10. Sample plots of the 3-dimensional joint angle curves estimated by employing all the five RFU methods, where "static" denotes the static RFU method [28], "linear" denotes the linearly interpolted RFU method [23], "joint axis" denotes the joint axis-based RFU method [29], "gravity" denotes the gravity acceleration based method [13] and "ours" denotes the RFU method proposed in this paper. Figure 11. RMSEs of each joint and each dimension averaged over trials and subjects, where "static" denotes the static RFU method [28], "linear" denotes the linearly interpolted RFU method [23], "joint axis" denotes the joint axis-based RFU method [29], "gravity" denotes the gravity acceleration based method [13] and "ours" denotes the RFU method proposed in this paper, * denotes the significant difference analyzed by the one-way ANOVA (α < 0.05). intra/extra rotation 6.8 ± 0.6 6.6 ± 0.5 2.3 ± 0.5 3.9 ± 0.4 2.0 ± 0.3 As shown in Table 2, the dispersion of correction quaternions of each RFU method range from 2.37 deg to 4.37 deg. The correction quaternion of linear interpolated RFU method presents the smallest dispersion, while the correction quaternion of joint axis-based RFU method presents the largest.

Discussion
This paper sought to firstly formulate the issue of reference frame unification and distinguish it with drift correction, and secondly to design, implement and evaluate a novel RFU method with comprehensive metrics based on the experimental analysis of the reference frames' deviations. We evaluate our method by estimating human lowerlimb joints from IMUs' measurements in the specially designed heterogeneous magnetic field and compare the performance with conventional RFU methods. The RMSEs of 3dimensional angle estimation that uses our RFU method are all under 3 deg, which is significantly better than the accuracy of using other methods. The standard deviations and dispersions of our method present similar repeatability with other methods, which as discussed in [23] suggests its promising application on clinical usage. These findings indicate that our proposed RFU method fits better to the time-varying characteristics of reference frames' deviations, thus measures and compensates for the deviations with state-of-the-art accuracy.
The accuracy of other RFU methods reported in our work is generally larger than that reported in literatures. The RMSEs on knee angles of the static, linearly interpolated and joint axis-based methods are all relatively higher than those reported in their original papers. This might result from two aspects. On one hand, as stated above, our experiments are conducted in the heterogeneous magnetic field, which as indicated in the study [30] significantly increase the estimation errors. On the other hand, studies of some literatures are performed using mechanical gimbals, as discussed in [23,34], the performance on the gimbal is different from that on human subjects, especially considering the fact that biological joints are neither ideal hinge nor screw joint. The RMSEs are also larger than those reported in [25] from which we employed the functional calibration postures. Although it is reported that the sensor-to-body alignment is magnetic distortion-free, the sensor-toreference frame orientation is sensitive to the magnetic distortion and other factors like movement-caused accelerations' corruption. In addition, the accuracy of the gravity acceleration-based method reported in our work is similar to that reported in its original work [13]. This is because our experiment is conducted for level walking, while the study [13] is conducted for skiing, which presents a larger range of motion and movementcaused accelerations.
The accuracy of our proposed RFU method is significantly better than that of others. The reason could be twofold. Firstly, as indicated by our analysis of the reference frames' deviations, the RFU method is supposed to be able to compensate for the non-linearly timevarying deviations. That is, the RFU method would better be pointwise. Compared with the static and linearly interpolated methods, our proposed method provides a pointwise compensation and the fusion coefficient is time-varied in each interval. Secondly, the metrics used to measure the deviations are more comprehensive than the metrics used in other methods. Compared with the joint axis-based and gravity acceleration-based methods that are also pointwise, our method's better performance demonstrates the efficiency of utilizing the more comprehensive metrics so that effects of the corruptions of each single metric can be decreased.
The repeatability and precision depicted by the dispersions and standard deviations present similarity with the data reported in literatures. The dispersions of the gravity acceleration-based method are slightly lower than that of its original study. It is because, as stated in [13], the original study performed on skiing results in a greater dispersion than normal lower-limb motions, like level walking. As for the dispersion of the joint axis-based RFU method, it is surprising that its dispersion is larger than those of static and linear methods, although it meets the larger precision reported in [29]. As discussed in [29], the hinge-joint constraint cannot always be satisfied, especially for the hip that is reported with greater similarity to the screw joint. The relatively larger dispersion and precision of our proposed method might result from the involvement of the gravity acceleration and joint axis.
The current study has several limitations. Firstly, performance is just evaluated on functional calibration-based sensor-to-body alignment and the accuracy is thus limited by the employed aligning method. There is an alternative way of aligning sensors to body segments without employing calibrations [27]. Further tests will be performed to evaluate the performance on such aligning methods. Secondly, the performance is also limited by the absolute orientation estimating performances of the Algorithm 1 we employed [22]. The reported performance only reflects the RFU method-caused accuracy improvement on this algorithm. Thirdly, performance is just evaluated on level walking. More motion patterns should be included for future work.

Conclusions
In this study, we clarify and formulate the issue of reference frame unification and propose a novel RFU method with a more comprehensive metric based on the analysis of the time-varying characteristics of reference frames' deviations. The analysis indicates the non-linear time-varying characteristics of deviations and the RFU's difference with drift correction, which further instructes the design of our method. The proposed method, which unifies the reference frames with a more comprehensive metric in a pointwise manner, is validated by the results. Specifically, the RMSEs of our method (2.7 deg, 1.4 deg, 1.1 deg for ankle, 2.6 deg, 1.5 deg, 1.3 deg for knee and 2.8 deg, 2.1 deg, 2.0 deg for hip) are significantly smaller than those of previously used RFU methods. These performances indicate that our method, even in the heterogeneous magnetic field, can achieve acceptable accuracy. The repeatability depicted by dispersion also demonstrates our method can be used across multiple usages with robust performances. This work would assist in further improving the accuracy of IMU-based joint angle estimation.

Institutional Review Board Statement:
The study was conducted according to the guidelines of the Declaration of Helsinki, and the protocol was approved by Chinese Ethics Committee of Registering Clinical Trials (ChiECRCT20200319).
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.