# Calibration of MEMS Triaxial Accelerometers Based on the Maximum Likelihood Estimation Method

Authors: Yifan Sun; Xiang Xu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

As a widely used inertial device, a MEMS triaxial accelerometer has zero-bias error, nonorthogonal error, and scale-factor error due to technical defects. Raw readings without calibration might seriously affect the accuracy of inertial navigation system. Therefore, it is necessary to conduct calibration processing before using a MEMS triaxial accelerometer. This paper presents a MEMS triaxial accelerometer calibration method based on the maximum likelihood estimation method. The error of the MEMS triaxial accelerometer comes into question, and the optimal estimation function is established. The calibration parameters are obtained by the Newton iteration method, which is more efficient and accurate. Compared with the least square method, which estimates the parameters of the suboptimal estimation function established under the condition of assuming that the mean of the random noise is zero, the parameters calibrated by the maximum likelihood estimation method are more accurate and stable. Moreover, the proposed method has low computation, which is more functional. Simulation and experimental results using the consumer low-cost MEMS triaxial accelerometer are presented to support the abovementioned superiorities of the maximum likelihood estimation method. The proposed method has the potential to be applied to other triaxial inertial sensors.


---

## Body

## 1. Introduction

Nowadays, with the gradual rise of microelectromechanical system (MEMS), low-precision inertial sensors, especially low-cost inertial sensors, have been widely used in many fields, such as motion tracking, attitude-controlling system, and unmanned aircraft systems. Of them, low-cost MEMS triaxial accelerometers are widely used in pedestrian navigation, vehicle navigation, unmanned aerial vehicle (UAV), carrier rocket, and other fields [1–5]. However, MEMS triaxial accelerometers have errors due to inaccuracy in manufacturing techniques, such as zero-bias error, scale-factor error, and nonorthogonal error. The existence of these errors will affect the accuracy of the raw readings of the MEMS triaxial accelerometer, thus affecting the accuracy of inertial navigation system. Therefore, for low-cost MEMS triaxial accelerometer sensors, calibration is a necessary step.Under laboratory conditions, the turntable is usually used to calibrate the MEMS triaxial accelerometer [6, 7]. However, MEMS triaxial accelerometer calibration is usually carried out without external auxiliary equipment in the field. In order to calibrate the MEMS triaxial accelerometer in the field test where it is inconvenient to use external auxiliary equipment such as a turntable, some methods have been proposed. In references [8–10], the 6-position method is used to collect the data of the MEMS triaxial accelerometer under the condition of static position and the error parameters are calibrated by the least square (LS) method. However, the scale-factor error and zero-bias error can be calibrated when the 6-position method is used, but the nonorthogonal error cannot be calibrated. Hence, a complete error model was constructed in references [11–13]. The LS method is used to estimate the error parameters, but these methods cannot obtain the optimal estimation because these methods assumed that the mean of the random noises is zero.Therefore, in order to improve the accuracy of calibration, the unbiased objective estimation function is constructed by the maximum likelihood (ML) estimation method in references [14–16]. In references [14, 15], the particle swarm optimization (PSO) algorithm is proposed to compute the objective function established by the ML estimation method, which is used to solve for the parameters. The initial value of iteration is not required, but the algorithm needs a longer time to obtain the calibrated results, so it is not suitable for practical application. In reference [16], the gradient descent method is proposed to estimate the error parameters, this process can be taken as the ML estimation problem. In each iteration, the error parameters are optimized by the gradient descent method. However, this method is relative to the initial rough values of the error parameters. If the initial rough values are not suitable, the estimation cannot get the global optimal solution. The step size and iteration times should also be taken into account when optimizing the objective function with the gradient descent method. If the sampling data are large, it is going to take a long time to compute in the gradient descent method.Due to the aforementioned drawbacks of the traditional methods, this paper proposes a new calibration method. Firstly, the total error model function is determined, and then the ML estimation method was used for the optimal estimation. In this paper, different from the abovementioned methods, the Newton iteration method is used to optimize the objective function established by the ML estimation method. The Newton iteration method has the advantages of fast convergence speed and global optimal solution. In practical application, the calibration of the MEMS triaxial accelerometer can be realized without external auxiliary equipment.


## 2. Error Analysis and Model Establishment of MEMS Triaxial Accelerometer

The raw readings of the MEMS triaxial accelerometer are affected by errors during manufacturing and installation; the errors are constructed by zero-bias error, scale-factor error and nonorthogonal error, and so on. The zero-bias error is the error caused by the MEMS triaxial accelerometer sensor outputs not being zero when the ideal outputs are zero, which is denoted byba in this paper. The scale-factor error is mainly caused by the different sensitivities of each sensor’s sensitive axis, which is represented by a diagonal matrix Cs in (1), where si (i = x, y, z) represents the scale-factor deviation of each axis [17]:(1)Cs=1+sx1+sy1+sz.The nonorthogonal error refers to the three sensitive axes that are not orthogonal to each other during the installation process of the MEMS triaxial accelerometer. A schematic diagram of the nonorthogonal error is given in Figure1, where the (Sx-Sy-Sz) frame is the orthogonal frame and the (X-Y-Z) frame is the nonorthogonal frame. The nonorthogonal error can be represented by the upper triangular matrix Cn in (2), which is calculated from the projection of the orthogonal frame to the nonorthogonal frame.(2)Cn=1sina−sinb0cosasinccosb00cosccosb.Figure 1
Schematic diagram of nonorthogonal error ((Sx-Sy-Sz) is the orthogonal frame; (X-Y-Z) is the nonorthogonal frame).The MEMS triaxial accelerometer is not only affected by its own error factors but also affected by the measurement noise [18]. In general, the raw readings of the MEMS triaxial accelerometer can be expressed as (3). In (3), ykn represents the local gravity acceleration, which is approximately considered unchanged in a certain area under local navigation frame (n-frame) [19], Cnb represents the rotation matrix from the navigation frame (n-frame) to the body frame (b-frame), and ε represents the Gaussian white noise with the variance of σ and the mean value of zero.(3)ykm=CSCnCnbykn+ba+ε.DefineTa=CSCn and ykb=Cnbykn. The error model (3) can be rewritten as (4), which is the MEMS triaxial accelerometer reading generation model with error:(4)ykm=Taykb+ba+ε.DefineRa=Ta−1, the error correction model of MEMS triaxial accelerometer under body frame (b-frame) is obtained as follows:(5)ykb=Raykm−ba−ε.


## 3. Calibration of Maximum Likelihood Estimation Method for MEMS Triaxial Accelerometer

### 3.1. Maximum Likelihood Estimation Model Construction

In this section, model (4) is used to construct the ML estimation model for parameters that need to be calibrated. In order to get the likelihood function Lθ, it is necessary to discuss the statistical properties of ykm. The random noise ε in model (4) is Gaussian white noise, which obeys the normal distribution. Then, the probability density function of ykm can be obtained as follows:(6)fykm=1σ2πexp−ykm−Taykb−ba22σ2.According to the theory of maximum likelihood estimation method, the likelihood function of a continuous-type random variable is(7)Lθ=∏k=1Nfykm,where fykm is the probability density function of sample ykm and θ is the parameter to be estimated, θ=Ta,ba,ykc. The idea of an ML estimate is to find the right θ such that the likelihood function is maximized:(8)maxθLθ=∏k=1Nfykm=1σ2πNexp∑k=1N−ykm−Taykb−ba22σ2.Finding the maximum value of likelihood function is equivalent to finding the minimum value of∑k=1Nykm−Taykb−ba2/2σ2; in order to facilitate the discussion, the constant term 2σ2 was ignored and then the ML problem is expressed as (9). Equation (10) can be obtained by expanding the constraint conditions in (9) using model (5). Normalizing the MEMS triaxial accelerometer calibrated data, the modular square of the data should be one. Now, due to the noise, the modular square of the calibrated MEMS triaxial accelerometer data is influenced.(9)minθ∑k=1Nek2=minθ∑k=1Nykm−Taykb−ba2,s.t.ykb2=1,k=1,…,(10)1=ykb2=Raykm−ba−ε2,=Raykm−ba2−2ykm−baTRaTRaε+εTRaTRaε,=ΔRaykm−ba2+γ,(11)μγ=ΔEγ=EtrRTRεεT=trERTRεεT,=trERTREεεT=trRTRσ2>0.According to (11), μγ is not equal to zero; it means that γ is a biased noise. If the noise γ is taken as the unbiased noise, like the traditional LS method, the estimated results are not optimal. In order to establish a total error estimation model to perform the optimal estimate, a Lagrange constraint term needs to be introduced to eliminate the influence of the biased noises γ. The minimization objective function is shown as (12), where λk denotes the Lagrange multiplier used to constrain the two paradigms of the calibration value ykc. The estimated value at this time is θ , θ=vecTTabaTy1bT⋯yNbTλ1⋯λNT, where vec· denotes the recomposition of a matrix into a column vector in the order of the columns, and the column vector does not include the three lower triangular elements of the matrix [20].(12)fθ=∑k=1Nykm−Taykb−ba2+λkykb2−1.


### 3.2. Newtonian Iteration Solution

The iteration formula for the Newtonian iteration method is(13)θi+1=θi−∇2fθi−1∇fθi,where θi+1 represents the (i + 1)-th iteration of parameter θ, θi represents the i-th iteration of parameter θ, ∇fθi and ∇2fθi represent the Jacobian matrix and Hessian matrix of the objective function in the i-th iteration process, respectively. The Jacobian matrix and Hessian matrix are described as follows:(14)∇f=∇fTaT∇fbaT∇fykbT︸k=1:N∇fλkT︸k=1:NT,∇2f=∇2fTaTa∇2fTaba∇2fTaykb⋯⋯∇2fTabaT∇2fbaba∇2fbaykb⋯⋯∇2fTaykbT∇2fbaykbT∇2fykbykb∇2fykbλk⋮⋮⋱⋱09×1T03×1T∇2fykbλkT0⋯⋮⋮⋱⋮,k=1:N.In the following formula, definedk=ykm−ba:(15)∇fTa=−2∑k=1Nykb⊗dk−Tykb,∇fba=−2∑k=1Ndk−Tykb,∇fykb=−2TaTdk−Tykb+2λkykb,∇fλk=ykb2−1,∇2fTaTa=2∑k=1NykbykbT⊗I,∇2fTaba=2∑k=1Nykb⊗I,∇2fTaykb=2ykb⊗IT−I⊗dk−Taykb,∇2fbaba=2NI,∇2fbaykb=2T,∇2fykbykb=2TTT+2λkI,∇2fykbλk=2ykb,where I denotes the 3×3 unit matrix and ⊗ denotes the Kronecker product. The abovementioned formulas are used in the Newton iteration method process. However, the Newton iteration method requires a good initial value to ensure that the iteration converges to a global minimum rather than a local minimum, so good initial estimate parameters are necessary for an accurate iteration result.


### 3.3. Initial Estimate

In this section, the ellipsoid fitting method of initial estimation is discussed. Since the initial parameters are required to ensure the accuracy before doing the ML estimate, it is important to give accurate initial parameters. Ignoring the noise in (10) and expanding 1=Raykm−ba2 give(16)ykmTEykm+2FTykm+G=0,where E=RaTRa, F=−RaTRaba, and G=baTRaTRaba−1. The output vectors of the MEMS triaxial accelerometer sensors lie on an ellipsoidal surface, and the coefficients of the ellipsoidal equation are a function of the error parameters, so an estimation of the parameters of calibration can be considered as an estimate of the parameters of the ellipsoid. The general equation of the ellipsoidal equation is(17)aXk2+bY2+cZ2+2fXY+2gXZ+2hYZ+2pX+2qY+2rZ+d=0.Bringing the vectorykm=XkYkZkT into (16) and then expanding and comparing with (17), we can get the following:(18)E=afgfbhghc,F=pqr,G=d.Substituting the raw readings into (17), the parameters of the ellipsoid can be obtained by the LS method, and then the initial parameters of calibration are obtained by Ra0=cholE and ba0=−E−1F. chol· is the Cholesky decomposition.


## 3.1. Maximum Likelihood Estimation Model Construction

In this section, model (4) is used to construct the ML estimation model for parameters that need to be calibrated. In order to get the likelihood function Lθ, it is necessary to discuss the statistical properties of ykm. The random noise ε in model (4) is Gaussian white noise, which obeys the normal distribution. Then, the probability density function of ykm can be obtained as follows:(6)fykm=1σ2πexp−ykm−Taykb−ba22σ2.According to the theory of maximum likelihood estimation method, the likelihood function of a continuous-type random variable is(7)Lθ=∏k=1Nfykm,where fykm is the probability density function of sample ykm and θ is the parameter to be estimated, θ=Ta,ba,ykc. The idea of an ML estimate is to find the right θ such that the likelihood function is maximized:(8)maxθLθ=∏k=1Nfykm=1σ2πNexp∑k=1N−ykm−Taykb−ba22σ2.Finding the maximum value of likelihood function is equivalent to finding the minimum value of∑k=1Nykm−Taykb−ba2/2σ2; in order to facilitate the discussion, the constant term 2σ2 was ignored and then the ML problem is expressed as (9). Equation (10) can be obtained by expanding the constraint conditions in (9) using model (5). Normalizing the MEMS triaxial accelerometer calibrated data, the modular square of the data should be one. Now, due to the noise, the modular square of the calibrated MEMS triaxial accelerometer data is influenced.(9)minθ∑k=1Nek2=minθ∑k=1Nykm−Taykb−ba2,s.t.ykb2=1,k=1,…,(10)1=ykb2=Raykm−ba−ε2,=Raykm−ba2−2ykm−baTRaTRaε+εTRaTRaε,=ΔRaykm−ba2+γ,(11)μγ=ΔEγ=EtrRTRεεT=trERTRεεT,=trERTREεεT=trRTRσ2>0.According to (11), μγ is not equal to zero; it means that γ is a biased noise. If the noise γ is taken as the unbiased noise, like the traditional LS method, the estimated results are not optimal. In order to establish a total error estimation model to perform the optimal estimate, a Lagrange constraint term needs to be introduced to eliminate the influence of the biased noises γ. The minimization objective function is shown as (12), where λk denotes the Lagrange multiplier used to constrain the two paradigms of the calibration value ykc. The estimated value at this time is θ , θ=vecTTabaTy1bT⋯yNbTλ1⋯λNT, where vec· denotes the recomposition of a matrix into a column vector in the order of the columns, and the column vector does not include the three lower triangular elements of the matrix [20].(12)fθ=∑k=1Nykm−Taykb−ba2+λkykb2−1.


## 3.2. Newtonian Iteration Solution

The iteration formula for the Newtonian iteration method is(13)θi+1=θi−∇2fθi−1∇fθi,where θi+1 represents the (i + 1)-th iteration of parameter θ, θi represents the i-th iteration of parameter θ, ∇fθi and ∇2fθi represent the Jacobian matrix and Hessian matrix of the objective function in the i-th iteration process, respectively. The Jacobian matrix and Hessian matrix are described as follows:(14)∇f=∇fTaT∇fbaT∇fykbT︸k=1:N∇fλkT︸k=1:NT,∇2f=∇2fTaTa∇2fTaba∇2fTaykb⋯⋯∇2fTabaT∇2fbaba∇2fbaykb⋯⋯∇2fTaykbT∇2fbaykbT∇2fykbykb∇2fykbλk⋮⋮⋱⋱09×1T03×1T∇2fykbλkT0⋯⋮⋮⋱⋮,k=1:N.In the following formula, definedk=ykm−ba:(15)∇fTa=−2∑k=1Nykb⊗dk−Tykb,∇fba=−2∑k=1Ndk−Tykb,∇fykb=−2TaTdk−Tykb+2λkykb,∇fλk=ykb2−1,∇2fTaTa=2∑k=1NykbykbT⊗I,∇2fTaba=2∑k=1Nykb⊗I,∇2fTaykb=2ykb⊗IT−I⊗dk−Taykb,∇2fbaba=2NI,∇2fbaykb=2T,∇2fykbykb=2TTT+2λkI,∇2fykbλk=2ykb,where I denotes the 3×3 unit matrix and ⊗ denotes the Kronecker product. The abovementioned formulas are used in the Newton iteration method process. However, the Newton iteration method requires a good initial value to ensure that the iteration converges to a global minimum rather than a local minimum, so good initial estimate parameters are necessary for an accurate iteration result.


## 3.3. Initial Estimate

In this section, the ellipsoid fitting method of initial estimation is discussed. Since the initial parameters are required to ensure the accuracy before doing the ML estimate, it is important to give accurate initial parameters. Ignoring the noise in (10) and expanding 1=Raykm−ba2 give(16)ykmTEykm+2FTykm+G=0,where E=RaTRa, F=−RaTRaba, and G=baTRaTRaba−1. The output vectors of the MEMS triaxial accelerometer sensors lie on an ellipsoidal surface, and the coefficients of the ellipsoidal equation are a function of the error parameters, so an estimation of the parameters of calibration can be considered as an estimate of the parameters of the ellipsoid. The general equation of the ellipsoidal equation is(17)aXk2+bY2+cZ2+2fXY+2gXZ+2hYZ+2pX+2qY+2rZ+d=0.Bringing the vectorykm=XkYkZkT into (16) and then expanding and comparing with (17), we can get the following:(18)E=afgfbhghc,F=pqr,G=d.Substituting the raw readings into (17), the parameters of the ellipsoid can be obtained by the LS method, and then the initial parameters of calibration are obtained by Ra0=cholE and ba0=−E−1F. chol· is the Cholesky decomposition.


## 4. Simulation and Experimental Analysis

### 4.1. Simulation Analysis

According to models (3) and (4), the parameters of the simulation are taken to be as (20). The attitude transformation matrix Cnb, which is converted from the navigation frame (n-frame) to the body frame (b-frame), is generated by Euler angles. The pitch angle, swing angle, and heading angle are taken to be 120, 120, and 80 degrees, respectively; to simulate the generation of real angle data, the data information was multiplied by random numbers to generate 300 sets of angle data, and then 300 sets of the attitude conversion matrix are converted. The random errors are set as Gaussian white noises with the mean of zero and a standard deviation of 0.1.(19)Ta=1.15000.0602−0.040101.09850.0575001.0978,ba=0.10.150.2.The ideal model of the vectors measured by the MEMS triaxial accelerometer after normalization should be a unit sphere, the center of the sphere is the coordinate origin, and the radius is 1. The model fitted from the measured and calibrated data is represented by the black ellipsoid and the data are represented by the red points in Figures2(a) and 2(b). It can be observed that the calibrated ellipsoid is a unit sphere with its center at the coordinate origin.Figure 2
(a) The ellipsoid fitted by measured data (black) and the data points (red). (b) The ellipsoid fitted by data calibrated by the ML estimation method (black) and the data points (red).
(a)(b)In Figure3, the MEMS triaxial accelerometer data calibrated by the LS method and the data calibrated by the ML estimation method are significantly improved compared with the uncalibrated MEMS triaxial accelerometer data. The modular square of the MEMS triaxial accelerometer calibrated by the ML estimation method fluctuates less and is closer to one than that of the LS method. Statistically, in Figure 3, the means of modular square are 1.172 (uncalibrated), 0.9945 (calibrated by the LS method), and 1 (calibrated by the ML method), respectively. The standard deviations of modular square are 0.18 (uncalibrated), 0.1052 (calibrated by the LS method), and 0.0403 (calibrated by the ML method), respectively. Both the LS method and the ML method have obvious effects on the calibration, but the results calibrated by the ML estimation method have less fluctuation in the vicinity of modular square of one than the results calibrated by the LS method. The parameters set in this paper and the parameters calibrated by the LS method and the ML estimation method are compared in Table 1. vecTa can be compared to observe the calibration effect of the LS method and ML estimation method for nonorthogonal error and scale-factor error, and ba can be compared to observe the calibration effect of the LS method and ML estimation method for zero-bias error. Compared with the simulation parameters, vecTa and ba calibrated by the LS method are reduced [−0.0827, 0.0066, 0.0298, 0.0096, 0.0165, and 0.1012] and [−0.131, 0.0062, and −0.0210], respectively. Similarly, compared with the simulation parameters, vecTa and ba calibrated by the ML estimation method are reduced [0.0077, 0.0073, −0.0007, −0.0013, −0.0063, and 0.0308] and [−0.0089, 0.0058, and −0.0079], respectively. By comparing with the parameters, it can be seen that both the LS method and ML estimation method can effectively calibrate the parameters. However, compared with the LS method, the parameters calibrated by the ML estimation method are closer to the simulation parameters, so it is more accurate.Figure 3
Comparison of uncalibrated (blue), LS calibrated (green), and ML estimation calibrated (pink) modular square.Table 1
Comparison of the parameters set in this article and the parameters calibrated by LS and ML estimation methods.
ParameterTabaSimulation parameter values1.15000.0602−0.040101.09850.0575001.09780.10.150.2Parameter values after calibration by the LS method1.06730.0668−0.030501.12830.0410001.19900.08690.15620.1790Parameter values after calibration by the ML method1.15770.0675−0.041401.09780.0512001.12860.09110.15580.1921In order to verify the stability of the ML estimation method under different random noise disturbances, 50 Monte Carlo experiments were conducted on the objective function values. The experiment results are shown in Figure4, which shows that the objective function values reach the convergence state after the third iteration. To further highlight the advantages of the ML estimation method over the LS method, the equations in (20) are defined as a metric and 50 Monte Carlo experiments are performed. In (20), Tacandbac, respectively, represent the parameters calibrated and Taandba, respectively, represent parameters set in this paper. In Figures 5 and 6, the variation of the objective function specified during the Monte Carlo experiment can be observed. The results of eT and eb, calculated by the parameters calibrated by the LS method, are both represented by the black solid lines. Also, the mean values of the results calibrated by the LS method are both represented by black dashed lines. Similarly, the results of eT and eb, calculated by the parameters calibrated by the ML method, are both represented by the red solid lines and the mean values of the results calibrated by the ML method are both represented by red dashed lines. In Figures 5 and 6, the parameters calibrated by the ML estimation method in the Monte Carlo experiment are closer to the parameters set in this paper than those calibrated by the LS method. In Figure 7, we can see that although the modular square of the data calibrated by the LS method is close to 1, the modular square of the data calibrated by the ML estimation method fluctuates slightly at 1. According to the analysis of Figures 5–7, it can be seen that the ML estimation method is more stable and accurate under the influence of noise than the LS method.(20)eT=16vecTa−Tac×100%,eb=13ba−bac×100%.Figure 4
50 Monte Carlo experiments of the objective function of the ML estimation method.Figure 5
50 Monte Carlo simulation experiments’eT change plot.Figure 6
50 Monte Carlo simulation experiments’eb change plot.Figure 7
Comparison of 50 Monte Carlo experiments of the mean of modular square calibrated by the LS method and the mean of modular square calibrated by the ML estimation method.


### 4.2. Experimental Analysis

This section will verify the proposed method through practical experiments. Figures8(a) and 8(b) show the experimental equipment and sensors of this paper, including the 3D printed frame, MPU9250 module, bubble level, and black desktop. The MPU9250 module is fixed to the 3D printed frame (the angle between the MPU9250 module and the 3D printed frame is not considered). The bubble level is used to observe the status of the black desktop. The sampling frequency of the MPU9250 module is 100 Hz, and the data are transmitted to the master computer through a serial port and then received and stored in the master computer. The 3D printed frame is placed on the black desktop. When the x-axis of the 3D frame is up, the 3D printed frame rotates horizontally by 90 degrees at regular intervals, so the static data of the accelerometer can be obtained when the x-axis is up. Similarly, static data can be obtained when the x-axis is down, when the y-axis is up and down, and when the z-axis is up and down. The sensor’s static data at 24 positions are recorded as data 1. The 3D printed frame is held in hand for slow movement, and the data obtained during the movement are recorded as data 2.Figure 8
(a, b) Experimental equipment and sensors.
(a)(b)Calibrating data 1 and data 2 by the LS method and ML estimation method, the modular square results of Figures9 and 10 are obtained, respectively. The statistical information can be obtained from Table 2. In Table 2, the mean of modular square of the two groups of sensors data calibrated by the LS method is very close to 1, while the mean of modular square calibrated by the ML method is 1, which is closer to the theoretical value 1, so it is more accurate. The standard deviation after calibration by the ML estimation method is significantly reduced compared with that calibrated by the LS method. In Figures 9 and 10, the pink curve fluctuates less than the green curve and fluctuates around 1.Figure 9
Comparison of modular square before and after calibrating measured data 1.Figure 10
Comparison of modular square before and after calibrating measured data 2.Table 2
Comparison of the standard deviation of modular square and mean of modular square of the uncalibrated data and the data calibrated by the LS method and ML method.
#Data 1#Data 2Uncalibrated mean of modular square1.1591.195Uncalibrated standard deviation of modular square0.24960.27LS method mean of modular square0.99960.9972LS method standard deviation of modular square0.028610.07423ML method mean of modular square11ML method standard deviation of modular square0.028530.07142The two sets of sensors data and calibration data are projected onto the ellipsoid, which is plotted in Figures11 and 12. By comparing Figures 11(a), 11(b), 12(a), and 12(b), after the calibration, the uncalibrated ellipsoid in black has been calibrated as a unit sphere, and the center of sphere is at the origin of the frame, which verifies the effectiveness of the calibration method proposed in this paper.Figure 11
(a) The ellipsoid fitted by data 1 (black) and the data 1 points (red). (b) The ellipsoid fitted by data 1 calibrated by the ML estimation method (black) and the points of data 1 calibrated by the ML estimation method (red).
(a)(b)Figure 12
(a) The ellipsoid fitted by data 2 (black) and the data 2 points (red). (b) The ellipsoid fitted by data 2 calibrated by the ML estimation method (black) and the points of data 2 calibrated by the ML estimation method (red).
(a)(b)


## 4.1. Simulation Analysis

According to models (3) and (4), the parameters of the simulation are taken to be as (20). The attitude transformation matrix Cnb, which is converted from the navigation frame (n-frame) to the body frame (b-frame), is generated by Euler angles. The pitch angle, swing angle, and heading angle are taken to be 120, 120, and 80 degrees, respectively; to simulate the generation of real angle data, the data information was multiplied by random numbers to generate 300 sets of angle data, and then 300 sets of the attitude conversion matrix are converted. The random errors are set as Gaussian white noises with the mean of zero and a standard deviation of 0.1.(19)Ta=1.15000.0602−0.040101.09850.0575001.0978,ba=0.10.150.2.The ideal model of the vectors measured by the MEMS triaxial accelerometer after normalization should be a unit sphere, the center of the sphere is the coordinate origin, and the radius is 1. The model fitted from the measured and calibrated data is represented by the black ellipsoid and the data are represented by the red points in Figures2(a) and 2(b). It can be observed that the calibrated ellipsoid is a unit sphere with its center at the coordinate origin.Figure 2
(a) The ellipsoid fitted by measured data (black) and the data points (red). (b) The ellipsoid fitted by data calibrated by the ML estimation method (black) and the data points (red).
(a)(b)In Figure3, the MEMS triaxial accelerometer data calibrated by the LS method and the data calibrated by the ML estimation method are significantly improved compared with the uncalibrated MEMS triaxial accelerometer data. The modular square of the MEMS triaxial accelerometer calibrated by the ML estimation method fluctuates less and is closer to one than that of the LS method. Statistically, in Figure 3, the means of modular square are 1.172 (uncalibrated), 0.9945 (calibrated by the LS method), and 1 (calibrated by the ML method), respectively. The standard deviations of modular square are 0.18 (uncalibrated), 0.1052 (calibrated by the LS method), and 0.0403 (calibrated by the ML method), respectively. Both the LS method and the ML method have obvious effects on the calibration, but the results calibrated by the ML estimation method have less fluctuation in the vicinity of modular square of one than the results calibrated by the LS method. The parameters set in this paper and the parameters calibrated by the LS method and the ML estimation method are compared in Table 1. vecTa can be compared to observe the calibration effect of the LS method and ML estimation method for nonorthogonal error and scale-factor error, and ba can be compared to observe the calibration effect of the LS method and ML estimation method for zero-bias error. Compared with the simulation parameters, vecTa and ba calibrated by the LS method are reduced [−0.0827, 0.0066, 0.0298, 0.0096, 0.0165, and 0.1012] and [−0.131, 0.0062, and −0.0210], respectively. Similarly, compared with the simulation parameters, vecTa and ba calibrated by the ML estimation method are reduced [0.0077, 0.0073, −0.0007, −0.0013, −0.0063, and 0.0308] and [−0.0089, 0.0058, and −0.0079], respectively. By comparing with the parameters, it can be seen that both the LS method and ML estimation method can effectively calibrate the parameters. However, compared with the LS method, the parameters calibrated by the ML estimation method are closer to the simulation parameters, so it is more accurate.Figure 3
Comparison of uncalibrated (blue), LS calibrated (green), and ML estimation calibrated (pink) modular square.Table 1
Comparison of the parameters set in this article and the parameters calibrated by LS and ML estimation methods.
ParameterTabaSimulation parameter values1.15000.0602−0.040101.09850.0575001.09780.10.150.2Parameter values after calibration by the LS method1.06730.0668−0.030501.12830.0410001.19900.08690.15620.1790Parameter values after calibration by the ML method1.15770.0675−0.041401.09780.0512001.12860.09110.15580.1921In order to verify the stability of the ML estimation method under different random noise disturbances, 50 Monte Carlo experiments were conducted on the objective function values. The experiment results are shown in Figure4, which shows that the objective function values reach the convergence state after the third iteration. To further highlight the advantages of the ML estimation method over the LS method, the equations in (20) are defined as a metric and 50 Monte Carlo experiments are performed. In (20), Tacandbac, respectively, represent the parameters calibrated and Taandba, respectively, represent parameters set in this paper. In Figures 5 and 6, the variation of the objective function specified during the Monte Carlo experiment can be observed. The results of eT and eb, calculated by the parameters calibrated by the LS method, are both represented by the black solid lines. Also, the mean values of the results calibrated by the LS method are both represented by black dashed lines. Similarly, the results of eT and eb, calculated by the parameters calibrated by the ML method, are both represented by the red solid lines and the mean values of the results calibrated by the ML method are both represented by red dashed lines. In Figures 5 and 6, the parameters calibrated by the ML estimation method in the Monte Carlo experiment are closer to the parameters set in this paper than those calibrated by the LS method. In Figure 7, we can see that although the modular square of the data calibrated by the LS method is close to 1, the modular square of the data calibrated by the ML estimation method fluctuates slightly at 1. According to the analysis of Figures 5–7, it can be seen that the ML estimation method is more stable and accurate under the influence of noise than the LS method.(20)eT=16vecTa−Tac×100%,eb=13ba−bac×100%.Figure 4
50 Monte Carlo experiments of the objective function of the ML estimation method.Figure 5
50 Monte Carlo simulation experiments’eT change plot.Figure 6
50 Monte Carlo simulation experiments’eb change plot.Figure 7
Comparison of 50 Monte Carlo experiments of the mean of modular square calibrated by the LS method and the mean of modular square calibrated by the ML estimation method.


## 4.2. Experimental Analysis

This section will verify the proposed method through practical experiments. Figures8(a) and 8(b) show the experimental equipment and sensors of this paper, including the 3D printed frame, MPU9250 module, bubble level, and black desktop. The MPU9250 module is fixed to the 3D printed frame (the angle between the MPU9250 module and the 3D printed frame is not considered). The bubble level is used to observe the status of the black desktop. The sampling frequency of the MPU9250 module is 100 Hz, and the data are transmitted to the master computer through a serial port and then received and stored in the master computer. The 3D printed frame is placed on the black desktop. When the x-axis of the 3D frame is up, the 3D printed frame rotates horizontally by 90 degrees at regular intervals, so the static data of the accelerometer can be obtained when the x-axis is up. Similarly, static data can be obtained when the x-axis is down, when the y-axis is up and down, and when the z-axis is up and down. The sensor’s static data at 24 positions are recorded as data 1. The 3D printed frame is held in hand for slow movement, and the data obtained during the movement are recorded as data 2.Figure 8
(a, b) Experimental equipment and sensors.
(a)(b)Calibrating data 1 and data 2 by the LS method and ML estimation method, the modular square results of Figures9 and 10 are obtained, respectively. The statistical information can be obtained from Table 2. In Table 2, the mean of modular square of the two groups of sensors data calibrated by the LS method is very close to 1, while the mean of modular square calibrated by the ML method is 1, which is closer to the theoretical value 1, so it is more accurate. The standard deviation after calibration by the ML estimation method is significantly reduced compared with that calibrated by the LS method. In Figures 9 and 10, the pink curve fluctuates less than the green curve and fluctuates around 1.Figure 9
Comparison of modular square before and after calibrating measured data 1.Figure 10
Comparison of modular square before and after calibrating measured data 2.Table 2
Comparison of the standard deviation of modular square and mean of modular square of the uncalibrated data and the data calibrated by the LS method and ML method.
#Data 1#Data 2Uncalibrated mean of modular square1.1591.195Uncalibrated standard deviation of modular square0.24960.27LS method mean of modular square0.99960.9972LS method standard deviation of modular square0.028610.07423ML method mean of modular square11ML method standard deviation of modular square0.028530.07142The two sets of sensors data and calibration data are projected onto the ellipsoid, which is plotted in Figures11 and 12. By comparing Figures 11(a), 11(b), 12(a), and 12(b), after the calibration, the uncalibrated ellipsoid in black has been calibrated as a unit sphere, and the center of sphere is at the origin of the frame, which verifies the effectiveness of the calibration method proposed in this paper.Figure 11
(a) The ellipsoid fitted by data 1 (black) and the data 1 points (red). (b) The ellipsoid fitted by data 1 calibrated by the ML estimation method (black) and the points of data 1 calibrated by the ML estimation method (red).
(a)(b)Figure 12
(a) The ellipsoid fitted by data 2 (black) and the data 2 points (red). (b) The ellipsoid fitted by data 2 calibrated by the ML estimation method (black) and the points of data 2 calibrated by the ML estimation method (red).
(a)(b)


## 5. Conclusions

In this paper, a complete error model of MEMS triaxial accelerometers is constructed. The influence of biased noises of the traditional objective function model, which is constructed by modulus calculation, is eliminated. To eliminate this influence, a calibration method based on the ML estimation method is proposed in this paper. The total error unbiased objective function is constructed, and the optimal estimation is obtained by the Newton iteration method. To make the Newton iteration method converge to the optimal value, an ellipsoid fitting method is adopted to get the initial parameters for the Newton iteration method. In the simulation, the results of the ellipsoid model and modular square verify the validity of the ML estimation method. Comparing using the Monte Carlo experiment, the results calibrated by the ML estimation method are better and more stable than the LS method. The experiment verifies the feasibility of the ML estimation method under different sets of data collected by the MPU9250 module. The method proposed in this paper is not only limited to calibrating the MEMS triaxial accelerometer but also has potential calibration effects for sensors of the same triaxial structure, such as magnetometers.



---

*Source: 4617365-2020-09-21.xml*